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

Ps.729

于 2024-06-22 09:58:55 发布

阅读量754

点赞数 18

文章标签： python matlab 开发语言

本文链接：https://blog.csdn.net/weixin_67304359/article/details/139870918

版权

💥💥💥💞💞💞欢迎来到本博客❤️❤️❤️💥💥💥
🎉作者研究：🏅🏅🏅主要研究方向是电力系统和智能算法、机器学习和深度学习。目前熟悉python网页爬虫、机器学习、群智能算法、深度学习的相关内容。希望将计算机和电网有效结合！⭐️⭐️⭐️

📋📋📋本文目录如下：⛳️⛳️⛳️

目录

0 概述

1 蒙特卡洛模拟方法介绍

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

3 流程图

4 运行结果

5 结果分析

6 参考文献

7 Matlab代码和Python代码

0 概述

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

1 蒙特卡洛模拟方法介绍

计算机模拟中的蒙特卡洛法也被称为随机抽样技术或统计检验方法，该方法最重要的特点是它是一种基于概率统计理论的方法。随着科学技术的发展和电子计算机的发明，蒙特卡洛法以其描述物理发展特点和物理实验过程的优点，在各个领域得到了广泛的应用。
蒙特卡洛法将目前所解决的问题当作是一种随机事件的概率，也可以将其看作是随机事件的期望值。按照某种实验的方式，通过某随机事件的出现频率来计算该事件的概率，或者是求得其数字特征，将实验所得的结论作为问题的解。采用蒙特卡洛法，在进行模拟计算时一般按照以下步骤进行,若所需解决的问题存在随机性特征，则就能将这个概率过程更加准确的描述和模拟，若问题并不具备随机性，反而是一种确定性问题，就必须认为建立相应的概率过程，而某些相关参量恰好就是该问题的解，整个过程就是将某问题的确定性改造为随机性。
在建立这种概率模型之后，可将不同的概率模型视为有不同类型的概率分布组成，由此就可获得己知概率分布的随机变量，如此就能成为可以进行蒙特卡罗方法模拟实验的主要方式，所以这种方法也是一种随机抽样。蒙特卡罗模拟的实现便是在已知概率分布所导致的随机数基础上进行的。
一般在建立这种概率模型后可以从中抽样，也就是实现模拟实验，在此之后还应当将其中某个随机变量确定下来，即所需解决问题的解，被称为无偏估计。通过获得多种估计量，例如方差等，就能够通过模拟方法获得问题的解。
对于本文的研究,依据不同用途电动汽车影响因素的分布函数和设定参数,采用蒙特卡洛法，对各用途电动汽车的日行驶里程、起始充电时间概率分布参数进行随机抽样，计算初始荷电状态和和充电所需时长，进而预测得出各类型的电动汽车充电负荷曲线，最后通过叠加各用途电动汽车的充电负荷曲线得出总的充电负荷曲线。

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

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

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

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

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

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

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

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

(2）依据不同用途电动汽车影响因素的分布函数和设定参数，如表1所示，采用蒙特卡洛模拟方法进行仿真，随机抽取日期类型、电动汽车起始充电时间和日行驶里程。
(3）计算车辆的初始荷电状态和充电所需时长。电动汽车动力电池的剩余电量直接关系到电动汽车的充电所需时间。对于电动汽车的充电所需时间，文章运用电动汽车动力电池的荷电状态( state of charge，SOC)来进行计算。假设电池消耗电量与行驶距离成正比， $d$ 为已行驶距离， $d_{m}$ 为电动模式下最大续航里程。此处假定每种不同用途的电动汽车每公里的耗电量相同,最后一次出行结束时的剩余电量由下式计算所得:【公式比较重要，纯手打】

$S O C_{1}=\left(S O C_{2}-\frac{d}{d_{m}}\right) \times 100 \%$
式中SOC1为完成充电时的电池荷电状态,SOC2为上一次完成充电时的电池荷电状态。

计算电动汽车充电所需时间T,可以通过将电池容量C、起始荷电状态SOC和充电功率Р求
得，具体公式如下:

$T=\frac{C(1-S O C)}{P}$
(4)计算某一种用途电动汽车在第i个充电负荷计算点时的总充电负荷。本文将每天计算为1440分钟，每15分钟计算一次充电负荷，共计96个充电负荷计算点，计算预测北京市各用途的电动汽车充电负荷，得到各用途电动汽车充电负荷。某一种用途电动汽车的充电负荷的预测计算方法如下:

$t_{1}=t_{2}+T$

其中 $t_{1}$ 表示第n 台某一种用途电动汽车结束充电的时刻， $t_{2}$ 表示第n台某种用途电动汽车开始充电的时刻，T表示第n台某种用途电动汽车充电所需时间。某一种用途电动汽车在第 i个充电负荷计算点时的总充电负荷 $P_{n,i}$ 可由以下方法计算得到
$P_{n, i}=\sum_{1}^{N} P_{t_{i}}$

$P_{t_{i}}$ 表示第n台电动汽车在第i个充电负荷计算点时的充电负荷，N表示某一种用途电动汽车的保有量。
(5)通过叠加各用途电动汽车的充电负荷得到总的电动汽车充电负荷。第i个充电负荷计算点的总电动汽车充电负荷的计算方式如公式(4-6)所示:

$P_{n, i}=\sum_{1}^{N} P_{t_{i}}P_{i}=\sum_{n c=1}^{N c} P_{n c, i}+\sum_{n t=1}^{N t} P_{n t, i}+\sum_{n b=1}^{N b} P_{n b, i}$

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

3 流程图

4 运行结果

部分代码：

%% 私家车负荷预测
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
%% ========结束==========

5 结果分析

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