基于蒙特卡洛法的规模化电动车有序充放电及负荷预测(Python和Matlab实现)

简介: 基于蒙特卡洛法的规模化电动车有序充放电及负荷预测(Python和Matlab实现)

0 概述

本文首先介绍蒙特卡洛模拟方法,再结合蒙特卡洛模拟方法、电动汽车充电负荷的影响因素和行驶数据及参数,确立基于蒙特卡洛法的规模化电动汽车充电负荷预测计算方法。根据北京市各用途的电动汽车未来保有量预测数据,对北京市各用途的电动汽车进行负荷预测,得到负荷曲线,并研究负荷曲线和对电网的影响结果。


1 蒙特卡洛模拟方法介绍

   计算机模拟中的蒙特卡洛法也被称为随机抽样技术或统计检验方法,该方法最重要的特点是它是一种基于概率统计理论的方法。随着科学技术的发展和电子计算机的发明,蒙特卡洛法以其描述物理发展特点和物理实验过程的优点,在各个领域得到了广泛的应用。

蒙特卡洛法将目前所解决的问题当作是一种随机事件的概率,也可以将其看作是随机事件的期望值。按照某种实验的方式,通过某随机事件的出现频率来计算该事件的概率,或者是求得其数字特征,将实验所得的结论作为问题的解。采用蒙特卡洛法,在进行模拟计算时一般按照以下步骤进行,若所需解决的问题存在随机性特征,则就能将这个概率过程更加准确的描述和模拟,若问题并不具备随机性,反而是一种确定性问题,就必须认为建立相应的概率过程,而某些相关参量恰好就是该问题的解,整个过程就是将某问题的确定性改造为随机性。

  在建立这种概率模型之后,可将不同的概率模型视为有不同类型的概率分布组成,由此就可获得己知概率分布的随机变量,如此就能成为可以进行蒙特卡罗方法模拟实验的主要方式,所以这种方法也是一种随机抽样。蒙特卡罗模拟的实现便是在已知概率分布所导致的随机数基础上进行的。

一般在建立这种概率模型后可以从中抽样,也就是实现模拟实验,在此之后还应当将其中某个随机变量确定下来,即所需解决问题的解,被称为无偏估计。通过获得多种估计量,例如方差等,就能够通过模拟方法获得问题的解。

  对于本文的研究,依据不同用途电动汽车影响因素的分布函数和设定参数,采用蒙特卡洛法,对各用途电动汽车的日行驶里程、起始充电时间概率分布参数进行随机抽样,计算初始荷电状态和和充电所需时长,进而预测得出各类型的电动汽车充电负荷曲线,最后通过叠加各用途电动汽车的充电负荷曲线得出总的充电负荷曲线。

2 规模化电动汽车充电负荷预测计算方法

根据对各用途电动汽车的充电影响因素进行研究得出了基本参数,包括充电时段起始充电容量分布和起始时间分布等,如表1所示。


该表给出了各用途电动汽车在建立充电负荷预测计算时需要的参数:


电动公交车一般进行常规充电,日行驶里程数和起始充电时间服从均匀分布;


电动出租车一般进行快速充电,日行驶里程数服从均匀分布,起始充电时间服从正态随机分布;


电动私家车进行常规充电和快速充电两种补给电能方式,日行驶里程数服从指数分布,起始充电时间服从正态随机分布。在常规充电和快速充电时,不同类型电动汽车的充电功率不尽相同,电动公交车的充电功率近似为电动私家车的五倍,电动出租车的充电功率近似为电动私家车的两倍。


经过原理及模型研究,假定各用途电动汽车均处于无序充电的状态,具体的计算方法如下:


(1)根据预测出的北京市各用途电动汽车的保有量,确定电动汽车的市场规模如表2所示。


(2)依据不同用途电动汽车影响因素的分布函数和设定参数,如表1所示,采用蒙特卡洛模拟方法进行仿真,随机抽取日期类型、电动汽车起始充电时间和日行驶里程。

(3)计算车辆的初始荷电状态和充电所需时长。电动汽车动力电池的剩余电量直接关系到电动汽车的充电所需时间。对于电动汽车的充电所需时间,文章运用电动汽车动力电池的荷电状态( state of charge,SOC)来进行计算。假设电池消耗电量与行驶距离成正比,为已行驶距离,为电动模式下最大续航里程。此处假定每种不同用途的电动汽车每公里的耗电量相同,最后一次出行结束时的剩余电量由下式计算所得:【公式比较重要,纯手打】


           

式中SOC1为完成充电时的电池荷电状态,SOC2为上一次完成充电时的电池荷电状态。

计算电动汽车充电所需时间T,可以通过将电池容量C、起始荷电状态SOC和充电功率Р求

得,具体公式如下:

               


(4)计算某一种用途电动汽车在第i个充电负荷计算点时的总充电负荷。本文将每天计算为1440分钟,每15分钟计算一次充电负荷,共计96个充电负荷计算点,计算预测北京市各用途的电动汽车充电负荷,得到各用途电动汽车充电负荷。某一种用途电动汽车的充电负荷的预测计算方法如下:

                 


其中表示第n 台某一种用途电动汽车结束充电的时刻,表示第n台某种用途电动汽车开始充电的时刻,T表示第n台某种用途电动汽车充电所需时间。某一种用途电动汽车在第 i个充电负荷计算点时的总充电负荷可由以下方法计算得到


       

表示第n台电动汽车在第i个充电负荷计算点时的充电负荷,N表示某一种用途电动汽车的保有量。

(5)通过叠加各用途电动汽车的充电负荷得到总的电动汽车充电负荷。第i个充电负荷计算点的总电动汽车充电负荷的计算方式如公式(4-6)所示:

   


式中,Nc,Nt,Nb分别表示在i时刻充电的电动私家车、电动出租车、电动公交车的数量;,,分别表示在i时刻电动私家车、电动出租车、电动公交车的充电负荷大小。

3 流程图

4 运行结果



5 Matlab代码和Python代码及文章

%% 私家车负荷预测
for n=1:1:N
    Pt1=zeros(1,1440);  %周一到周日
    for i=1:1:N1
        Ts=unifrnd( 1,7,1,1);
        Ts=round(Ts);
        %% ============工作日情况============
        if(Ts>=1&&Ts<=5) %Mon---Fri
            ds=exprnd(50.00,1,1);  %指数分布X~E (0.020);1/0.02=50[期望E(x)]
            if(ds>=180)
                ds=180;
            end
            if(ds<=0)
                ds=0;
            end
            soc1=1-(ds./dm(1)); %(4-1)
            Ts1=unifrnd(1,1440,1,1);
            %% ============工作日情况快充===========
            if(Ts1>=420&&Ts1<=1080) %时间7:00--18:00
                ts2=unifrnd(420,1080,1,1);  %均匀分布:[420,1080]
                ts2=round(ts2); %对充电时长取整
                ts3=( 1-soc1).*40*60/Ps(2);    %T
                ts3=round(ts3); %对充电时长取整
                k=ts2+ts3;  %(4-3)t1=t2+T
                if k>1440   %如果充电结束时间大于24时
                    k=k-1440;
                    Pt1(1:k)=Pt1(1:k)+Ps(2); %# 清晨充电
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(2); %晚上24时前充电
                else      %如果充电结束时间没有过24时,无需在清晨充电
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(2);  % 充电时段加上充电功率
                end
                %% ============工作日情况慢充===========
            else   %时间18:00--7:00
                ts2=normrnd( 1150,100,1,1);  %N (1150,100^2)
                ts2=round(ts2);
                if(ts2<=1)
                    ts2=1;
                end
                ts1=(1-soc1).*40*60/Ps(1);
                ts1=round(ts1);
                k=ts1+ts2;
                if k>1440
                    k=k-1440;
                    Pt1(1:k)=Pt1(1:k)+Ps(1);
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(1);
                else
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(1);
                end
            end
        end
        %% =======非工作日情况=========
        if(Ts>=6&&Ts<=7)
            ds=exprnd(55.55,1,1);  %X~E(0.018);1/0.018=55.55
            if (ds>=180)
                ds=180;
            end
            if (ds<=0)
                ds=0;
            end
            soc1=1-ds./dm(1);
            Ts1=unifrnd(1,1440,1,1);
            %% ======非工作日情况快充=====
            if(Ts1>=420&&Ts1<=1080)  %时间7:00--18:00
                ts2=unifrnd(420,1080,1,1);  %均匀分布:[420,1080]
                ts2=round(ts2);
                ts3=(1-soc1).*40*60/Ps(2);
                ts3=round(ts3);
                k=ts2+ts3;
                if k>1440
                    k=k-1440;
                    Pt1(1:k)=Ptl(1:k)+Ps(2);
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(2);
                else
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(2);
                end
                %% ======非工作日情况慢充=========
            else   %时间18:00--7:00
                ts2=normrnd(1150,150,1,1);
                ts2=round(ts2);
                if (ts2<=1)
                    ts2=l;
                end
                ts1=(1-soc1).*40*60/Ps(1);
                ts1=round(ts1);
                k=ts1+ts2;
                if k>1440
                    k=k-1440;
                    Pt1(1:k)=Pt1(1:k)+Ps(1);
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(1);
                else
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(1);
                end
            end
        end
        Pt1(1:1440)=Pt1(1:1440);
    end
    %% 公交车负荷预测
    Pt2=zeros(1,1440);
    for i=1:1:N2
        dg=unifrnd(150,200,1,1);  %日行驶里程数
        %% ========上下限限制======
        if(dg>=200)
            dg=200;
        end
        if (dg<-0)
            dg=0;
        end
        soc2=1-(dg./dm(2));
        tg1=(1-soc1).*300*60/Pg;
        tg1=round(tg1);
        tg2=unifrnd(1,1440,1,1);
        tg2=round(tg2);
        kg=tg1+tg2;
        if (kg>1440)
            kg=kg-1440;
            Pt2(1:kg)=Pt2(1:kg)+Pg;
            Pt2(tg2:1440)=Pt2(tg2:1440)+Pg;
        else
            Pt2(tg2:kg)=Pt2(tg2:kg)+Pg;
        end
        Pt2(1:1440)=Pt2(1:1440);
    end
    %% ========结束==========
%% 私家车负荷预测
for n=1:1:N
    Pt1=zeros(1,1440);  %周一到周日
    for i=1:1:N1
        Ts=unifrnd( 1,7,1,1);
        Ts=round(Ts);
        %% ============工作日情况============
        if(Ts>=1&&Ts<=5) %Mon---Fri
            ds=exprnd(50.00,1,1);  %指数分布X~E (0.020);1/0.02=50[期望E(x)]
            if(ds>=180)
                ds=180;
            end
            if(ds<=0)
                ds=0;
            end
            soc1=1-(ds./dm(1)); %(4-1)
            Ts1=unifrnd(1,1440,1,1);
            %% ============工作日情况快充===========
            if(Ts1>=420&&Ts1<=1080) %时间7:00--18:00
                ts2=unifrnd(420,1080,1,1);  %均匀分布:[420,1080]
                ts2=round(ts2); %对充电时长取整
                ts3=( 1-soc1).*40*60/Ps(2);    %T
                ts3=round(ts3); %对充电时长取整
                k=ts2+ts3;  %(4-3)t1=t2+T
                if k>1440   %如果充电结束时间大于24时
                    k=k-1440;
                    Pt1(1:k)=Pt1(1:k)+Ps(2); %# 清晨充电
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(2); %晚上24时前充电
                else      %如果充电结束时间没有过24时,无需在清晨充电
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(2);  % 充电时段加上充电功率
                end
                %% ============工作日情况慢充===========
            else   %时间18:00--7:00
                ts2=normrnd( 1150,100,1,1);  %N (1150,100^2)
                ts2=round(ts2);
                if(ts2<=1)
                    ts2=1;
                end
                ts1=(1-soc1).*40*60/Ps(1);
                ts1=round(ts1);
                k=ts1+ts2;
                if k>1440
                    k=k-1440;
                    Pt1(1:k)=Pt1(1:k)+Ps(1);
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(1);
                else
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(1);
                end
            end
        end
        %% =======非工作日情况=========
        if(Ts>=6&&Ts<=7)
            ds=exprnd(55.55,1,1);  %X~E(0.018);1/0.018=55.55
            if (ds>=180)
                ds=180;
            end
            if (ds<=0)
                ds=0;
            end
            soc1=1-ds./dm(1);
            Ts1=unifrnd(1,1440,1,1);
            %% ======非工作日情况快充=====
            if(Ts1>=420&&Ts1<=1080)  %时间7:00--18:00
                ts2=unifrnd(420,1080,1,1);  %均匀分布:[420,1080]
                ts2=round(ts2);
                ts3=(1-soc1).*40*60/Ps(2);
                ts3=round(ts3);
                k=ts2+ts3;
                if k>1440
                    k=k-1440;
                    Pt1(1:k)=Ptl(1:k)+Ps(2);
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(2);
                else
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(2);
                end
                %% ======非工作日情况慢充=========
            else   %时间18:00--7:00
                ts2=normrnd(1150,150,1,1);
                ts2=round(ts2);
                if (ts2<=1)
                    ts2=l;
                end
                ts1=(1-soc1).*40*60/Ps(1);
                ts1=round(ts1);
                k=ts1+ts2;
                if k>1440
                    k=k-1440;
                    Pt1(1:k)=Pt1(1:k)+Ps(1);
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(1);
                else
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(1);
                end
            end
        end
        Pt1(1:1440)=Pt1(1:1440);
    end
    %% 公交车负荷预测
    Pt2=zeros(1,1440);
    for i=1:1:N2
        dg=unifrnd(150,200,1,1);  %日行驶里程数
        %% ========上下限限制======
        if(dg>=200)
            dg=200;
        end
        if (dg<-0)
            dg=0;
        end
        soc2=1-(dg./dm(2));
        tg1=(1-soc1).*300*60/Pg;
        tg1=round(tg1);
        tg2=unifrnd(1,1440,1,1);
        tg2=round(tg2);
        kg=tg1+tg2;
        if (kg>1440)
            kg=kg-1440;
            Pt2(1:kg)=Pt2(1:kg)+Pg;
            Pt2(tg2:1440)=Pt2(tg2:1440)+Pg;
        else
            Pt2(tg2:kg)=Pt2(tg2:kg)+Pg;
        end
        Pt2(1:1440)=Pt2(1:1440);
    end
    %% ========结束==========
%% 私家车负荷预测
for n=1:1:N
    Pt1=zeros(1,1440);  %周一到周日
    for i=1:1:N1
        Ts=unifrnd( 1,7,1,1);
        Ts=round(Ts);
        %% ============工作日情况============
        if(Ts>=1&&Ts<=5) %Mon---Fri
            ds=exprnd(50.00,1,1);  %指数分布X~E (0.020);1/0.02=50[期望E(x)]
            if(ds>=180)
                ds=180;
            end
            if(ds<=0)
                ds=0;
            end
            soc1=1-(ds./dm(1)); %(4-1)
            Ts1=unifrnd(1,1440,1,1);
            %% ============工作日情况快充===========
            if(Ts1>=420&&Ts1<=1080) %时间7:00--18:00
                ts2=unifrnd(420,1080,1,1);  %均匀分布:[420,1080]
                ts2=round(ts2); %对充电时长取整
                ts3=( 1-soc1).*40*60/Ps(2);    %T
                ts3=round(ts3); %对充电时长取整
                k=ts2+ts3;  %(4-3)t1=t2+T
                if k>1440   %如果充电结束时间大于24时
                    k=k-1440;
                    Pt1(1:k)=Pt1(1:k)+Ps(2); %# 清晨充电
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(2); %晚上24时前充电
                else      %如果充电结束时间没有过24时,无需在清晨充电
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(2);  % 充电时段加上充电功率
                end
                %% ============工作日情况慢充===========
            else   %时间18:00--7:00
                ts2=normrnd( 1150,100,1,1);  %N (1150,100^2)
                ts2=round(ts2);
                if(ts2<=1)
                    ts2=1;
                end
                ts1=(1-soc1).*40*60/Ps(1);
                ts1=round(ts1);
                k=ts1+ts2;
                if k>1440
                    k=k-1440;
                    Pt1(1:k)=Pt1(1:k)+Ps(1);
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(1);
                else
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(1);
                end
            end
        end
        %% =======非工作日情况=========
        if(Ts>=6&&Ts<=7)
            ds=exprnd(55.55,1,1);  %X~E(0.018);1/0.018=55.55
            if (ds>=180)
                ds=180;
            end
            if (ds<=0)
                ds=0;
            end
            soc1=1-ds./dm(1);
            Ts1=unifrnd(1,1440,1,1);
            %% ======非工作日情况快充=====
            if(Ts1>=420&&Ts1<=1080)  %时间7:00--18:00
                ts2=unifrnd(420,1080,1,1);  %均匀分布:[420,1080]
                ts2=round(ts2);
                ts3=(1-soc1).*40*60/Ps(2);
                ts3=round(ts3);
                k=ts2+ts3;
                if k>1440
                    k=k-1440;
                    Pt1(1:k)=Ptl(1:k)+Ps(2);
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(2);
                else
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(2);
                end
                %% ======非工作日情况慢充=========
            else   %时间18:00--7:00
                ts2=normrnd(1150,150,1,1);
                ts2=round(ts2);
                if (ts2<=1)
                    ts2=l;
                end
                ts1=(1-soc1).*40*60/Ps(1);
                ts1=round(ts1);
                k=ts1+ts2;
                if k>1440
                    k=k-1440;
                    Pt1(1:k)=Pt1(1:k)+Ps(1);
                    Pt1(ts2:1440)=Pt1(ts2:1440)+Ps(1);
                else
                    Pt1(ts2:k)=Pt1(ts2:k)+Ps(1);
                end
            end
        end
        Pt1(1:1440)=Pt1(1:1440);
    end
    %% 公交车负荷预测
    Pt2=zeros(1,1440);
    for i=1:1:N2
        dg=unifrnd(150,200,1,1);  %日行驶里程数
        %% ========上下限限制======
        if(dg>=200)
            dg=200;
        end
        if (dg<-0)
            dg=0;
        end
        soc2=1-(dg./dm(2));
        tg1=(1-soc1).*300*60/Pg;
        tg1=round(tg1);
        tg2=unifrnd(1,1440,1,1);
        tg2=round(tg2);
        kg=tg1+tg2;
        if (kg>1440)
            kg=kg-1440;
            Pt2(1:kg)=Pt2(1:kg)+Pg;
            Pt2(tg2:1440)=Pt2(tg2:1440)+Pg;
        else
            Pt2(tg2:kg)=Pt2(tg2:kg)+Pg;
        end
        Pt2(1:1440)=Pt2(1:1440);
    end
    %% ========结束==========


完整代码:基于蒙特卡洛法的规模化电动车有序充放电及负荷预测

6 结果分析

研究分析北京市各用途电动汽车充电负荷曲线,我们发现北京市总充电负荷曲线类似于北京市电动私家车的充电负荷曲线,正是由于电动私家车在总电动汽车保有量中占比最大,对总负荷的影响最大。而电动公交车的充电起始时间较为均匀,因此电动公交车的充电负荷曲线分布较为平均,曲线动荡主要来自不同站点的出发时间和到站时间的不同。电动出租车的充电负荷主要集中在电动出租车的常规充电时间节点上,因此电动出租车的充电负荷曲线波动也处于常规充电时间的节点上。


7 参考文献

[1]郑晶晶,闫志杰,李伟,王伟,彭飞云,董海鹰.基于蒙特卡洛法的电动汽车无序充电对电网的影响分析[J].电气传动自动化,2019,41(05):1-5.


[2]蒋林洳,万伟江,丁霄寅,李涛永,张元星,张晶.一种基于直接蒙特卡洛法的电动汽车充电负荷模型[J].供用电,2018,35(04):20-25+13.DOI:10.19421/j.cnki.1006-6357.2018.04.004.


[3]李冰轩. 基于蒙特卡洛法的规模化电动汽车充电负荷预测研究[D].宁夏大学,2018.


[4]陈鹏,孟庆海,赵彦锦.基于蒙特卡洛法的电动汽车充电负荷计算[J].电气工程学报,2016,11(11):40-46.

目录
打赏
0
0
0
0
77
分享
相关文章
基于CVX凸优化的电动汽车充放电调度matlab仿真
本程序基于CVX凸优化实现电动汽车充放电调度,通过全局和局部优化求解,展示了不同情况下的负载曲线。程序在MATLAB 2022a上运行,有效平抑电网负荷峰值,提高电网稳定性。
|
4月前
|
基于python-django的matlab护照识别网站系统
基于python-django的matlab护照识别网站系统
30 0
|
6月前
|
【Python】实现MATLAB中计算两个矩形相交面积的rectint函数
Python中实现MATLAB中rectint函数的方法,该函数用于计算两个矩形相交区域的面积,并通过定义Rectangle类和calc_area函数展示了如何计算两个矩形的交集面积。
88 1
基于蒙特卡洛的电力系统可靠性分析matlab仿真,对比EDNS和LOLP
电力系统可靠性评估研究,聚焦于LOLP(电力不足概率)和EDNS(期望缺供电量)的模拟分析。使用MATLAB2022a进行基于蒙特卡洛的仿真,模拟单线及多线故障,分析连锁效应。程序中通过随机断开线路,计算潮流,判断越限并用PSO优化。结果显示,LOLP和EDNS增加时,故障概率降低,但小概率大影响事件概率上升。以IEEE24-RTS系统为案例,考虑元件失效的马尔科夫过程,不考虑3个及以上元件失效情况,因为可能导致系统大规模崩溃。仿真步骤包括随机线路断开、故障分析和稳定性评估,涉及信息节点概率计算、潮流计算及优化决策。
探索Python编程:从基础到高级
在这篇文章中,我们将一起深入探索Python编程的世界。无论你是初学者还是有经验的程序员,都可以从中获得新的知识和技能。我们将从Python的基础语法开始,然后逐步过渡到更复杂的主题,如面向对象编程、异常处理和模块使用。最后,我们将通过一些实际的代码示例,来展示如何应用这些知识解决实际问题。让我们一起开启Python编程的旅程吧!
Python编程入门:从零基础到实战应用
本文是一篇面向初学者的Python编程教程,旨在帮助读者从零开始学习Python编程语言。文章首先介绍了Python的基本概念和特点,然后通过一个简单的例子展示了如何编写Python代码。接下来,文章详细介绍了Python的数据类型、变量、运算符、控制结构、函数等基本语法知识。最后,文章通过一个实战项目——制作一个简单的计算器程序,帮助读者巩固所学知识并提高编程技能。
[oeasy]python053_学编程为什么从hello_world_开始
视频介绍了“Hello World”程序的由来及其在编程中的重要性。从贝尔实验室诞生的Unix系统和C语言说起,讲述了“Hello World”作为经典示例的起源和流传过程。文章还探讨了C语言对其他编程语言的影响,以及它在系统编程中的地位。最后总结了“Hello World”、print、小括号和双引号等编程概念的来源。
126 80

推荐镜像

更多
AI助理

你好,我是AI助理

可以解答问题、推荐解决方案等