模拟大规模电动车充电行为（Matlab实现）

1 模拟大规模充电汽车充电行为

2 Matlab部分代码

3 Matlab代码实现 

1 模拟大规模充电汽车充电行为

电动汽车EV(Electric Vehicle)具有清洁环保、高效节能的优点,不仅能缓解化石能源危机,而且能够有效地减少温室气体的排放。2015年10月，国务院发布加快EV充电基础设施建设的指导意见,指出到2020年充电基础设施能满足500万辆EV充电需求,预计未来几年我国EV的保有量将大幅增长。然而,规模化EV的无序充电会加大电网负荷的峰谷差,并对电力系统的规划、配电网的电能质量和经济运行以及稳定性带来显著的影响，反之.对EV的充电行为进行有序优化控制,充分发挥EV作为分布式储能元件的优势,能够实现削峰填谷、平抑可再生能源出力波动的功能,并为电网提供调峰、调频等辅助服务。

下图为大规模EV分散接入配电网的场景示意图。EV的管理框架分为配电网代理商、本地代理商和EV 3个层次。配电网代理商分布在高中压HVIMV(High VoltageMedium Voltage)变电层中,负责配电网的安全稳定控制以及EV的有序协调控制;本地代理商分布在中低压MVLV(Mediuim VoltageLow Voltage)变电层中,负责区域EV充电负荷的管控。充电站连接在配电变压器下,配电变压器下除了EV负荷EVL(EV Load)外,还有常规的居民负荷RL(Residential Load)。

EV进入充电站后,充电智能终端可以获取电池的总容量、荷电状态SoC(State Of Charge)等信息,并通过以太网或专用无线网络提交给本地代理商;用户设置取车时刻以及充电预期荷电状态,充电结束后按实际充电电量向本地代理商支付费用。本地代理商将EV的状态信息和电量需求信息汇集后上传给配电网代理商。配电网代理商采取分时段实时滚动优化的控制策略,当有新的EV接入电网时,更新EV的充电需求信息,执行优化,控制算法,并将充电计划分区下达给各个本地代理商,由本地代理商执行对管控区域内EV的充电控制。

2 Matlab部分代码实现

function [will_charge, best_prior, best_idx] = ...CalChargingPrior(is_necessary, cur_soc, arrMCost, ...pindex_start, pindex_end)global TOU_EPriceglobal solutionsglobal battery_features% debug by feiif (pindex_end - pindex_start + 1) < battery_features.fcharge_periodserrorend% 优先权重系数if solutions.start_charging == 4    %如果是组合策略，用上所有系数W = solutions.w_coeff;% 剩余SOC/最低SOCu_s = 1.5 - cur_soc / battery_features.lowest_soc;elseW = zeros(1, 4);W(solutions.start_charging) = 1;W(4) = 1;end% 是否必须充电if is_necessaryu_s = 10000;    %infelseu_s = 0;endnPrior = length(arrMCost);% 从某个时刻开始充电的平均电价：元/（瓦＊时段）%arrMCost = [];% TOC平均电价：元/（瓦＊时段）nMTOC = TOU_EPrice.mprice;% 充电需用时段nTch = single( battery_features.fcharge_periods );% 剩余空闲时段nTidle = pindex_end - pindex_start + 1;arrTidle = single(nTidle) : -1 : nTch;% 随机优先 [0 1]nRndPr = round(rand);arrRndPr = single( nRndPr * rand(1, nPrior) );% 计算不同时段的充电优先程度arrPrior = W(1)*arrRndPr + ...W(2)*(1.5 - arrMCost ./ nMTOC) + ...W(3)*(1.5 - nTch ./ arrTidle) + ...u_s + W(4);% 最优充电级别与时段[p, i] = max(arrPrior);best_prior = p;best_idx = pindex_start + i - 1;% 是否建议充电will_charge = (best_prior > 0.5);
end

% 根据不同策略，计算起始充电时段
function [start_cperiod, end_cperiod] = CalChargingPeriods(is_necessary, cur_soc, ...pindex_start, pindex_end)%global solutionsglobal battery_features%global mc_params%global TOU_EPrice% 初始参数start_cperiod = 0;end_cperiod = 0;earliest = pindex_start;latest = pindex_end - battery_features.fcharge_periods + 1;% 计算从某一时段起充电所需平均电价arrMCost = CalChargePrices(earliest:latest);% 根据策略，计算每一时段充电优先级[will_charge, best_prior, best_idx] = ...CalChargingPrior(is_necessary, cur_soc, arrMCost, pindex_start, pindex_end);% 计算策略下的最佳充电起止时间if will_chargestart_cperiod = best_idx;end_cperiod = start_cperiod + battery_features.fcharge_periods - 1;endend

% m_w_load: 一天中电网各时段负荷平均值，periods_per_day X mean_load
% m_w_charged: 一天中EV各时段充电平均值，periods_per_day X mean_charge
% m_real_trv：一天中EV各时段出行次数平均值，periods_per_day X mean_real_travel
% m_plan_trv：一天中EV各时段计划出行次数平均值，periods_per_day X mean_plan_travel
function [Y, Rmp, Rsc, Rt] = CalIndices(m_w_load, m_w_charged, ...m_real_trv, m_plan_trv, is_print)global TOU_EPrice%     W1 + W2 + W3 = 1.0
%     W1 * 0.275 = W2 * 0.440 = W3 * 0.975 = Si
%
%     Mi = [0.275 0.440 0.975]
%     Wi = [0.5455    0.3409    0.1538]
%     Si = [0.15 0.15 0.15]% 指标统计范围 及 权重系数RAGmp = [0.05 0.25];    RAGsave=[0.15 0.75];   RAGtrip=[0.95 1.0];Wi = [0.4 0.3 0.3];% 谷峰比 ＝ 充电谷值／充电峰值
%     nWValley = min(m_w_load);
%     nWPeak = max(m_w_load);
%     Rpv = nWValley / nWPeak;% 均峰比 ＝ 充电均值／充电峰值nWMeanLoad = mean(m_w_load);nWPeakLoad = max(m_w_load);Rmp = nWMeanLoad / nWPeakLoad;% 用户节省电费率 ＝ 1 - 实际电费/最高电费nRealCost = sum(m_w_charged .* TOU_EPrice.day_prices);nMaxCost = sum(m_w_charged .* max(TOU_EPrice.day_prices));Rsc = 1 - nRealCost / nMaxCost;% 顺利出行率 ＝ 实际里程/计划里程nRealTrv = sum(m_real_trv);nPlanTrv = sum(m_plan_trv);Rt = nRealTrv / nPlanTrv;%     Rmp = 0.218;
%     Rsc = 0.193;
%     Rt = 0.972;% 综合指标：均峰比，节能率，出行率Y = Wi(1) * (Rmp - min(RAGmp)) / (max(RAGmp) - min(RAGmp)) + ...Wi(2) * (Rsc - min(RAGsave)) / (max(RAGsave) - min(RAGsave)) + ...Wi(3) * (Rt - min(RAGtrip)) / (max(RAGtrip) - min(RAGtrip));if isnan(Y)Y = 0.0;end% 输出结果if is_printfprintf('----电网负荷：峰值= %d 瓦；均值= %d 瓦；均峰比= %f----\n', nWPeakLoad, nWMeanLoad, Rmp);fprintf('----居民用电费用：实际电费= %f 元；最高电费= %d 元；节省率= %f----\n', nRealCost, nMaxCost, Rsc);fprintf('----顺利出行情况：实际出行里程= %f 公里；计划出行里程= %d 公里；出行率= %f----\n', nRealTrv, nPlanTrv, Rt);fprintf('----综合指标：Y= %f----\n', Y);end
end

function [all_indices, mday_indices] = EVPowerLoad(varargin)global solutionsglobal mc_params%% 初始化变量[nTerms, nEVs] = Initial(varargin);%% 模拟 n天 m辆电动汽车if ~mc_params.only_show_resultsfor t=1:nTerms
%           for ev_id=1:total_EVall_indices = Simulate(t);
%             endendend%% 结果统计: 峰谷差，出行受影响概率，充电价格，加入噪声出行等result = TempResult('load');[grid, ev_behaviour, mday_indices] = StatisticPowerLoad(result);
end%% initial parameters
function [nTerms, nEVs] = Initial(params)%% 定义常量global mc_params            % MC 参数global battery_features     % 电池特性global pile_power           % 电桩特性global solutions            % 策略global behaviours           % 一天用户行为（EVs X Day_Periods）%global all_behaviours     % 所有统计时间的用户行为（EVs X All_Periods）global t_periods            % 所有统计时段global g2v_features           % 充电特性global TOU_EPrice         % 分时电价global PDF_Travel           % 服从一定概率分布的出行行为%% 指定参数运行if size(params, 2) > 0nTerms = params{1};nDaysPT = params{2};nEVs = params{3};chargingStrategy = params{4};chargingWCoeff = params{5};chargingMinSOC = params{6};isOutput = params{7};% 用于存储中间结果，而不用保存成文件saveResult = false;showResultOnly = false;else  % default values%     clear allclcclose all% 自定义参数值nTerms = 4;nDaysPT = 25;nEVs = 240000;             %EV数量: 2800 5000 62500 240000chargingStrategy = 4;chargingWCoeff = [0.1512,0.7384,0.1105,0.0196];%ones(1, 4) * 0.25;chargingMinSOC = 0.5686;%0.2;isOutput = true;saveResult = true;         %default:true;showResultOnly = true;      %default:false;seed=15;randn('state',seed);rand('state',seed);end%% 初始化常量mc_params.output = isOutput;                            %是否输出结果mc_params.cur_day = 1;                                 %当前天数%mc_params.cur_ev = 1;                                   %当前汽车IDmc_params.periods_per_day = 96;                           %一天划分为若干时段mc_params.total_days = nDaysPT;                               %每期模拟天数mc_params.all_days = nTerms * mc_params.total_days;       %总共模拟天数mc_params.mins_per_period = (24 * 60) / ...mc_params.periods_per_day;  %每个时段的时长mc_params.total_periods = mc_params.periods_per_day * ...(mc_params.total_days + 1);       %时段总数mc_params.total_EVs = nEVs;mc_params.cperiod_start_id = 0;mc_params.cperiod_end_id = 0;mc_params.cperiods = [];mc_params.eday_soc = ones(nEVs, 1, 'single');          %初始所有EV的SOC为1.0mc_params.save_result = saveResult;mc_params.memory_result = [];                           %用于存储中间结果，而不用保存成文件mc_params.only_show_results = showResultOnly;           %只输出统计结果battery_features.capacity = 100;                        %电池容量battery_features.voltage = 230;                         %电池电压battery_features.power = battery_features.capacity * ...battery_features.voltage;      %电池总能量battery_features.efficiency = 0.3;                      %充电效率battery_features.fcharge_duration = 5;                  %充满时间battery_features.full_soc = 0.9;                        %soc > 此状态，即认为电池充满battery_features.lowest_soc = chargingMinSOC;           %soc < 此状态，必须充电battery_features.power_consume_per_km = 0.125e3;        %耗能% 满充所需分钟数 及 时段数battery_features.fcharge_minutes = battery_features.fcharge_duration * 60;battery_features.fcharge_periods = ceil( battery_features.fcharge_minutes / mc_params.mins_per_period );pile_power = 15e3;                                      %充/放电功率solutions.start_charging = chargingStrategy;            %起始充电策略：1.随机充电；2.电价引导充电；3.停车即充电；4.组合策略solutions.w_coeff = chargingWCoeff;                     %各种充电影响因素的权重solutions.enable_discharge = false;                     %允许V2G的放电模式% 电动汽车类型%EV_type = {'bus', 'taxis', 'official_car', 'private_car'};% 1 X (模拟天数 * 每天时段数)t_periods = CreatePeriods();% 计算每时段电网耗用能量 (Wh)g2v_features.grid_w_consumed = mc_params.mins_per_period / 60 * pile_power;% 计算每时段EV所充能量 (Wh)g2v_features.ev_w_charged = g2v_features.grid_w_consumed * battery_features.efficiency;% 加载分时电价表TOU_EPrice = CreateEPriceList();% 清除结果数据if ~mc_params.only_show_resultsClearResults();end% 加载概率密度函数load './data/fitness_travel.mat';PDF_Travel = pdf_travel;
end%% simulate one day
function all_indices = Simulate(t)global mc_paramsglobal behaviours%global t_periodsmc_params.cur_day = t;%mc_params.cur_ev = ev_id;%w_day = mod( (d-1), 7 );%% clear all behavioursClearAllBehaviours();%fprintf('－－－－－模拟第 %d 天－－－－－\n', mc_params.cur_day);%% make the driving & charging plan for all daysMakeDaysPlan();%% update all EVsUpdateEVs();%% formating data[grid, ev_behaviour, y_indices] = StatDataFormat();%% combine data of everydayall_indices = CombineData(grid, ev_behaviour, y_indices);%% update paramsUpdateMCParams();
endfunction all_periods = CreatePeriods()global mc_params% construct the period arrayall_periods = int32( 0 : (mc_params.total_periods - 1) ) * ...mc_params.mins_per_period;
endfunction ClearAllBehaviours()global behaviours% init array for EV behavioursbehaviours = [];behaviours.v_is_driving = logical([]);behaviours.v_plan_driving = logical([]);behaviours.v_driving_km_pp = single([]);behaviours.v_driving_cost_power = single([]);behaviours.v_able_charge = logical([]);behaviours.v_is_charging = logical([]);behaviours.v_ev_w_charged = single([]);behaviours.v_ev_w_discharged = single([]);behaviours.grid_power_load = single([]);behaviours.soc = single([]);
endfunction ClearResults()global mc_paramsglobal solutionsndays = mc_params.all_days;% init & save temporary dataresult.grid.power_load = zeros(ndays, mc_params.periods_per_day, 'single');result.grid.power_charged = zeros(ndays, mc_params.periods_per_day, 'single');result.grid.power_V2G = zeros(ndays, mc_params.periods_per_day, 'single');result.ev_behaviour.trv_num_per_period = zeros(ndays, mc_params.periods_per_day, 'int32');result.ev_behaviour.plan_trv_num_per_period = zeros(ndays, mc_params.periods_per_day, 'int32');result.ev_behaviour.cha_num_per_period = zeros(ndays, mc_params.periods_per_day, 'int32');result.ev_behaviour.dischar_num_per_period = zeros(ndays, mc_params.periods_per_day, 'int32');result.y_indices = zeros(ndays, 1, 'single');TempResult('save', result);clear 'result';
endfunction all_indices = CombineData(grid, ev_behaviour, y_indices)global behaviours%global all_behavioursglobal solutionsglobal mc_paramsndays = mc_params.all_days;cur_day = ((mc_params.cur_day - 1) * mc_params.total_days + 1) : ...mc_params.cur_day * mc_params.total_days;% 加载临时数据result = TempResult('load');% 合并每天的电网负荷及EV行为result.grid.power_load(cur_day, :) = grid.power_load;result.grid.power_charged(cur_day, :) = grid.power_charged;result.grid.power_V2G(cur_day, :) = grid.power_V2G;result.ev_behaviour.trv_num_per_period(cur_day, :) = ev_behaviour.trv_num_per_period;result.ev_behaviour.plan_trv_num_per_period(cur_day, :) = ev_behaviour.plan_trv_num_per_period;result.ev_behaviour.cha_num_per_period(cur_day, :) = ev_behaviour.cha_num_per_period;result.ev_behaviour.dischar_num_per_period(cur_day, :) = ev_behaviour.dischar_num_per_period;result.y_indices(cur_day, :) = y_indices;% save temporary dataTempResult('save', result);% clear data from memoryall_indices = result.y_indices;clear 'result';
endfunction UpdateMCParams()global mc_paramsglobal behaviours% clear mc pramsmc_params.total_periods = mc_params.periods_per_day * ...(mc_params.total_days + 1);       %时段总数mc_params.cperiod_start_id = 0;mc_params.cperiod_end_id = 0;mc_params.cperiods = [];% record soc for all EVs when one day is endingmc_params.eday_soc = behaviours.soc(:, end);
end

function [is_necessary, soc_start] = IsNecessaryCharge(ev_id, pindex_start, pindex_end)%global behavioursglobal battery_features% soc低于0.2时必须充电soc_start = GetPreviousSOC(ev_id, pindex_start);is_necessary = (soc_start < battery_features.lowest_soc);
end

3 Matlab代码实现

完整Matlab代码实现：