1 主要内容

该程序参考《基于动态交通信息的电动汽车充电负荷时空分布预测》和《基于动态交通信息的电动汽车充电需求预测模型及其对配网的影响分析》文献模型，考虑私家车、出租车和共用车三类交通工具特性和移动负荷特性，实现了基于动态交通信息的电动汽车充电负荷时空分布预测。将负荷预测情况和33节点配电网络进行结合，形成交通网-配电网交互模型，采用牛拉法进行潮流计算。程序采用matlab编写，注释清晰，方便学习！

交通网-配电网交互模型

路网模型参考下图：

很多文献都会用到路网结合模型，这也是电动车研究领域的一个热点方向，通过路网模型更能精确定位电能需求点以及对配电网络的影响。

动态交通路网模型

动态交通网络模型采用图论法进行建模，拓扑结构示意图如下：

具体的数学模型表达式为：

式中：G 为交通路网；V 表示图 G 所有节点的集合；E 表示图 G 所有路段的集合；W 为路段权值的集合，即道路路阻；K 表示划分的时间段集合，即将全天划分为 m个时间段。

交通网络生成代码如下：

%%动态交通路网模型：采用图论分析方法建模
%G为交通路网（图）；V表示图G所有节点的集合，本文设置32个配电网节点；E表示所有路段的集合；W为路段权重的集合，即道路路阻；K表示划分的时间段集合，即将全天划分为m个时间段。
%路段权值W表示道路出行代价，可采用路段长度、通行速度、行程时间等权值进行量化研究
%在城市路网中，路口交叉节点多设置信号灯进行管控，车辆行驶既受到路段阻抗影响，又在交叉节点产生时间延误。
%因此，城市道路路阻可表示为W(t)=CV(t)+RV(t);CV(t)表示节点阻抗模型，RV(t)表示路段阻抗模型
%依据城市交通状况划分标准，饱和度 S 评价指标：畅通(0＜S≤0.6)、缓行(0.6＜S≤0.8)、拥挤(0.8＜S≤1.0)以及严重拥堵(1.0＜S≤2.0)。
%道路交叉口和路段通行能力不同，可以得到不同饱和度对应的路段阻抗和节点阻抗模型
%RV(t)路段阻抗模型，通过饱和度S=Q/C，（Q为路段交通流量，C为通行能力，这里是随机生成S所以不用管Q与C的问题）；t0为零流量行程时间；alpha beta为阻抗影响因子，来计算。
%CV(t)节点阻抗模型，通过信号周期c，绿信比lamda，路段车辆到达率q，来计算
%针对RV(t)路段阻抗和CV(t)节点阻抗模型中，饱和度S是唯一变量，其余为道路规划固定参数，因此，将RV(t)和CV(t)合并可得到道路路阻阻抗模型W(t)=CV(t)+RV(t)。

%% alpha beta 阻抗影响因子 t0-零流量行程时间，计算RV(t)路段阻抗模型
alpha=1.3;
beta=1.2;
t0=10;

%% 信号周期c，绿信比lamda，路段车辆到达率q，计算CV(t)节点阻抗模型
c=30;
lamda=0.7;
q=0.8;

%% 2【路网结构】 邻接矩阵(此时建立的仅是结构，而为给每条边赋以路段权值W)
LJ=zeros(32,32);
L=[1 2;2 3;3 9;1 4;2 6;3 7;4 5;5 6;6 7;7 8;8 9;4 10;5 11;6 12;7 13;10 11;11 12;12 13;13 14;8 21;9 26;10 15;11 16;12 17;12 18;14 20;15 16;16 17;17 18;18 19;19 20;20 21;15 32;16  27;17 22;18 23;19 24;20 24;21 25;32 27;27 22;22 23;23 24;24 25;25 26;23 28;24 29;25 30;26 31;27 28;28 29;29 30;30 31];
% 矩阵序号法写邻接矩阵
LJ(L(:,2)*32-32+L(:,1))=1; 
LJ(L(:,1)*32-32+L(:,2))=1;
%
figure(1)
G=graph(LJ,'upper');%根据带权邻接矩阵生成无向图
% plot(G);
% title('无向图')
plot(G,'EdgeLabel',G.Edges.Weight)
title('标定权重的无向图')

%% 3【阻抗模型权重W】(为给每条边赋以路段权值W)
% 1） 随机生成S权重   
% 2） 或者 输入shortestpath函数求路径（这一项操作没有便要做，不知道作用是为何）
S=zeros(32,32);
S=2*rand(32,32);
for i=1:32for j=1:32if LJ(i,j)==1if S(i,j)<=1 Rv(i,j)=t0*(1+alpha*S(i,j)^beta);elseRv(i,j)=t0*(1+alpha*(2-S(i,j))^beta);endif S(i,j)<=0.6 Cv(i,j)=0.9*(c*(1-lamda)^2/2/(1-lamda*S(i,j))+S(i,j)^2/2/q/(1-S(i,j)));elseCv(i,j)=c*(1-lamda)^2/2/(2-lamda*S(i,j))+1.5*(S(i,j)-0.6)*S(i,j)/(2-S(i,j));endW(i,j)=Rv(i,j)+Cv(i,j);elseW(i,j)=0*1e3;endend
end
W(find(W>50))=50;
figure(2)
bar3(W)
%
G=graph(W,'upper');
% plot(G);
% title('无向图')
plot(G,'EdgeLabel',G.Edges.Weight)
title('标定权重的无向图')

%[path,distance]=shortestpath(G,1,6)      %这行代码在这块没用，是多余的
%% 4【路程模型权重L】（作为计算耗电量中的一个参数）
RL=10+10*rand(32,32);
RL=RL.*LJ;

2 部分代码

load dataTRN.mat;
load dataEV.mat;
load dataWq.mat;
G=graph(W,'upper');                %生成路-网耦合图G，W是每条边的权重值，W为32*32阶矩阵
Charge=zeros(size(Mcar,1),5);      %生成一个size(Mcar,1)=1000*5的零矩阵，其中size(Mcar,1)=1000。

TP_carnumber=zeros(size(W,1),3);   %生成一个size(W,1)=32*3的零矩阵，W是每条边的权重值，W为32*32阶矩阵,TP_carnumber(x,1)、TP_carnumber(x,2)、TP_carnumber(x,3)应该是每条边上/每个交通节点私家车1，公交车2，私家车3的数量
Pcharge=zeros(24,3);               %生成一个24*3的零矩阵，******PCharge代表24小时内私家车1，公交车2，私家车3
Pchar_slow=12;                     %设置慢充充电功率为12KW
Pchar_fast=48;                     %设置慢充充电功率为48KW
Pntcharge=zeros(24,32);            %生成一个24*32的零矩阵，PntCharge代表24个小时内32个配网节点的供电负荷=电动汽车耗电量
%% MC
for icar=1:size(Mcar,1)             %EV的编号，从1-1000，共1000台EV，size(Mcar,1)=1000，一个for icar=1:size(Mcar,1) 循环代表把第一台EV的%1-9个特征变量（车辆编号icar、EV种类Icar_kind、EV初始位置Birthland、%目的地Destination、初始时刻tBirth1、返程时刻tBirth2、ev容量Cbat、初始soc SOC0、速度Vcar）全部采集一遍%1000个for icar=1:size(Mcar,1) 循环代表把1000台EV的%1-9个特征变量（车辆编号icar、EV种类Icar_kind、EV初始位置Birthland、%目的地Destination、初始时刻tBirth1、返程时刻tBirth2、ev容量Cbat、初始soc SOC0、速度Vcar）全部采集一遍icarIcar_kind=Mcar(icar,2);%EV种类；%size(Mcar,2)=1000Birthland=Mcar(icar,3);%EV初始位置Destination=Mcar(icar,4);%目的地tBirth1=Mcar(icar,5);%初始时刻tBirth2=Mcar(icar,6);%返程时刻Cbat=Mcar(icar,7);%ev容量SOC0=Mcar(icar,8);%初始socVcar=Mcar(icar,9);%速度TP_carnumber(Birthland,Icar_kind)=TP_carnumber(Birthland,Icar_kind)+1;%第Birthland条边/第Birthland个交通节点下，对应的第Icar_kind种车型加1，此行代码是计数的作用if   Mcar(icar,3)==Mcar(icar,4)                                       %如果EV初始位置Birthland=car(icar,3)==EV目的地Destination=Mcar(icar,4)continue                                                           %end                                                                   %[path,distance]=shortestpath(G,Birthland,Destination);                %如果EV初始位置Birthland=car(icar,3)   不等于  EV目的地Destination=Mcar(icar,4)%则输出EV初始位置Birthland=car(icar,3)与EV目的地Destination=Mcar(icar,4)最短路径path与最短距离distance%% C1 私家车if  Icar_kind==1%考虑环境温度和速度 耗电量dE=1.5*distance/Vcar+(0.21-1e-3*Vcar+1.531/Vcar)*distance;       %由最短行驶距离distance，电动汽车行驶实时车速Vcar=Mcar(icar,9)，计算出EV初始位置到目的地的耗电量SOC1=SOC0-dE/0.9/Cbat;                                           %由初始soc即SOC0=Mcar(icar,8) 与 EV初始位置到目的地的耗电量 与 充电效率0.9 与 ev容量Cbat=Mcar(icar,7)SOC2=SOC1-dE/0.9/Cbat;                                           %计算出电动汽车剩余电量SOC2tdest1=tBirth1+distance/Vcar;                             %到达目的地时间（tdest1）=初始时刻+行驶时间if  SOC2<0                                                %判断私家车是否充电的条件% 触发慢充/计及排队时间Tmmc=Wq(ceil(tdest1));                                %ceil为向上取整T80=tdest1+(0.8-SOC1)*Cbat/Pchar_slow+Tmmc;           %充电到0.8EV容量时间（T80）=到达目的地时间+慢充到0.8EV容量时间+Tmmcif T800.8EV容量时间＜到达第二个目的地的时间Charge(icar,2)=1;                                  %Icar_kind=Mcar(icar,2);EV种类；所以Charge(icar,2)赋值为私家车Icar_kind==1Charge(icar,3)=floor(tdest1);                      %EV初始位置Mcar(icar,3)赋值为向下取整tdest1，即32个交通节点中的一个Charge(icar,4)=ceil(T80);                          %目的地Mcar(icar,4)赋值为向上取整T80，即32个交通节点中的一个Charge(icar,5)=ceil(Pchar_slow);                   %初始时刻Mcar(icar,5)赋值为向上取整Pchar_slowPntcharge(Charge(icar,3):Charge(icar,4),Destination)=Pntcharge(Charge(icar,3):Charge(icar,4),Destination)+Charge(icar,5);%列表切片，Pntcharge(Charge(icar,3):Charge(icar,4),Destination)即为取Charge(icar,3):Charge(icar,4)所在行与Destination所在列的交叉元素elseCharge(icar,2)=1;                                  %Icar_kind=Mcar(icar,2);EV种类；所以Charge(icar,2)赋值为私家车Icar_kind==1Charge(icar,3)=floor(tdest1);                      %EV初始位置Mcar(icar,3)赋值为向下取整tdest1Charge(icar,4)=ceil(tBirth2);                      %目的地Mcar(icar,4)赋值为向上取整tBirth2，因为此时T80＞tBirth2成立Charge(icar,5)=ceil(Pchar_slow);                   %初始时刻Mcar(icar,5)赋值为向上取整Pchar_slowPntcharge(Charge(icar,3):Charge(icar,4),Destination)=Pntcharge(Charge(icar,3):Charge(icar,4),Destination)+Charge(icar,5);%列表切片，Pntcharge(Charge(icar,3):Charge(icar,4),Destination)即为取Charge(icar,3):Charge(icar,4)所在行与Destination所在列的交叉元素