✅作者简介：热爱科研的Matlab仿真开发者，擅长毕业设计辅导、数学建模、数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎 往期回顾关注个人主页：Matlab科研工作室

👇 关注我领取海量matlab电子书和数学建模资料

🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信。

🔥 内容介绍

一、光储充换电站发展背景

1. 电动汽车产业兴起

   ：随着全球对环境保护和可持续发展的重视，电动汽车作为一种清洁能源交通工具，其保有量迅速增长。这推动了充换电基础设施的大规模建设，光储充换电站作为一种综合性的充换电设施，集光伏发电、储能、充电与电池更换功能于一体，不仅能为电动汽车提供电能补给，还能有效整合分布式能源，提升能源利用效率。

2. 能源转型需求

   ：传统能源面临枯竭以及环境压力增大的问题，促使能源结构向可再生能源转型。太阳能作为一种丰富的清洁能源，在光储充换电站中得到广泛应用。通过光伏发电，为站内设备及电动汽车充电提供绿色能源，减少对传统电网的依赖，助力能源可持续发展。

二、用户充电负荷与分时电价互动的重要性

1. 削峰填谷

   ：电动汽车充电负荷具有随机性和集中性特点，若大量电动汽车在同一时段充电，会给电网带来巨大的负荷压力，可能导致电网峰谷差进一步增大。分时电价策略通过在不同时段设置不同的电价，引导用户在电价低谷时段充电，高峰时段减少充电行为，从而实现削峰填谷，平衡电网负荷，提高电网运行的稳定性和经济性。

2. 用户成本与体验优化

   ：对于用户而言，合理的分时电价政策使他们能够根据自身需求和电价变化，灵活安排充电时间，降低充电成本。同时，这种互动模式也能避免因充电过于集中而导致的等待时间过长等问题，提升用户的充电体验。

三、考虑互动的光储充换电站优化模型构建要素

1. 光伏发电模块

   ：光储充换电站中的光伏发电系统通过太阳能电池板将太阳能转化为电能。其发电功率受到光照强度、温度等因素影响，具有间歇性和波动性。在优化模型中，需要考虑实时的光照数据以及光伏发电系统的效率特性，预测光伏发电量，为站内提供绿色能源，并合理安排其在不同时段的使用。

2. 储能系统

   ：储能系统在光储充换电站中起到能量缓冲和调节作用。当光伏发电量过剩或电价较低时，储能系统充电储存能量；当光伏发电不足、用电高峰或电价较高时，储能系统放电，为充电设备供电。优化模型需考虑储能系统的充放电功率限制、充放电效率、容量衰减等因素，合理规划储能系统的充放电策略，以平衡站内能量供需，降低运营成本。

3. 充电与换电服务

   ：充换电站为电动汽车提供充电和换电两种服务模式。不同的充电方式（快充、慢充）功率不同，对电网负荷和用户等待时间的影响也不同。换电服务则涉及电池储备、更换效率等问题。在优化模型中，要综合考虑用户对充电和换电服务的需求、站内充电设备和电池储备情况，合理安排充电和换电资源，提高服务质量和运营效率。

4. 用户充电负荷响应

   ：用户充电行为受分时电价影响，会根据电价高低调整充电时间和充电量。优化模型需建立用户充电负荷响应模型，通过分析用户的充电习惯、用电需求弹性等因素，预测不同分时电价策略下用户的充电负荷变化，从而实现与分时电价的有效互动。

四、优化模型的目标与约束

1. 目标函数

   ：

   • 经济效益最大化

     ：光储充换电站运营商的主要目标之一是实现经济效益最大化。这包括光伏发电收益、向电网售电收益、充电和换电服务收益，同时减去购买电网电能成本、储能系统建设与运维成本、充电设备维护成本等运营成本。

   • 电网负荷平衡

     ：通过合理安排充电、换电、光伏发电及储能系统的运行，减小电网峰谷差，优化电网负荷曲线，降低对电网的冲击，提高电网运行稳定性，间接产生经济效益和社会效益。

2. 约束条件

   ：

   • 功率平衡约束

     ：站内光伏发电功率、储能系统充放电功率、充电设备消耗功率以及与电网交互功率之间需满足功率平衡关系，确保站内能量供需平衡。

   • 设备运行约束

     ：光伏发电系统、储能系统、充电设备等都有各自的技术参数和运行限制，如光伏发电功率上限、储能系统充放电功率限制、充电设备的额定功率等，模型需保证各设备在安全可行的范围内运行。

   • 电池管理约束

     ：对于换电服务，要考虑电池的储备数量、电池寿命、充电状态等因素，确保有足够的可用电池满足用户换电需求，同时合理管理电池充电过程，延长电池使用寿命。

⛳️ 运行结果

📣 部分代码

%% -- values

% parameters

mu          = 5;

sigma       = sqrt(0.5);

% time series

time_nod    = 500;

time        = linspace(mu,mu + 10*sigma,time_nod);

% nu values

nu_nod      = 11;

n           = 1;

nu          = linspace(0.01 + n - 1,0.9 + n - 1,nu_nod);

% nu            = [0.05:.1:0.9,0.95:0.01:0.99];%[0.05:0.1:0.95];

% nu_nod      = numel(nu);

%% -- calculations

%

% traditional solution

df1dt1      = @(t) (exp(-(mu - t).^2/(2*sigma^2)).*(2*mu - 2*t))/(2*sigma^3*sqrt(2*pi));

df2dt2      = @(t) -(exp(-(mu - t).^2/(2*sigma^2)).*(- mu^2 + 2*mu*t + sigma^2 - t.^2))/sigma^5/sqrt(2*pi);

df3dt3      = @(t) -(exp(-(mu - t).^2/(2*sigma^2)).*(mu - t).*(- mu^2 + 2*mu*t + 3*sigma^2 - t.^2))/sigma^7/sqrt(2*pi);

df4dt4      = @(t) (exp(-(mu - t).^2/(2*sigma^2)).*(mu^4 - 4*mu^3*t - 6*mu^2*sigma^2 + 6*mu^2*t.^2 + 12*mu*sigma^2*t - 4*mu*t.^3 + 3*sigma^4 - 6*sigma^2*t.^2 + t.^4))/sigma^9/sqrt(2*pi);

df5dt5      = @(t) (exp(-(mu - t).^2/(2*sigma^2)).*(mu - t).*(mu^4 - 4*mu^3*t - 10*mu^2*sigma^2 + 6*mu^2*t.^2 + 20*mu*sigma^2*t - 4*mu*t.^3 + 15*sigma^4 - 10*sigma^2*t.^2 + t.^4))/sigma^11/sqrt(2*pi);

df6dt6      = @(t) -(exp(-(mu - t).^2/(2*sigma^2)).*(- mu^6 + 6*mu^5*t + 15*mu^4*sigma^2 - 15*mu^4*t.^2 - 60*mu^3*sigma^2*t + 20*mu^3*t.^3 - 45*mu^2*sigma^4 + 90*mu^2*sigma^2*t.^2 - 15*mu^2*t.^4 + 90*mu*sigma^4*t - 60*mu*sigma^2*t.^3 + 6*mu*t.^5 + 15*sigma^6 - 45*sigma^4*t.^2 + 15*sigma^2*t.^4 - t.^6))/sigma^13/sqrt(2*pi);

%  

% -(t - mu).*exp(-(t - mu).^2/2/sigma^2)/sqrt(2*pi)/sigma^3;

dG1         = df1dt1(time);

dG2         = df2dt2(time);

dG3         = df3dt3(time);

dG4         = df4dt4(time);

dG5         = df5dt5(time);

dG6         = df6dt6(time);

% Gauss-Legendre Approach

% fractional solution

DG          = nan(nu_nod,time_nod);

DGn         = nan(nu_nod,time_nod);

for iinu  = 1 : nu_nod

    nux         = nu(iinu) - floor(nu(iinu));

    gamma_nu    = nux/(1 - nux);

    switch floor(nu(iinu))

        case 0

            DG(iinu,:)    = cfgaussder(time,nu(iinu),mu,sigma);

            [DGn(iinu,:),Par(iinu),new_time(iinu,:)] = cfgaussder_distr (time,nu(iinu),mu,sigma);

        case 1

            DG(iinu,:)  = (gamma_nu + 1)*dG1 - gamma_nu*cfgaussder(time,nux,mu,sigma);

        case 2

            DG(iinu,:)  = (gamma_nu + 1)*((-gamma_nu)*dG1 + dG2) + (- gamma_nu)^2*cfgaussder(time,nux,mu,sigma);

        case 3

            DG(iinu,:)  = (gamma_nu + 1)*((-gamma_nu)^2*dG1 + (-gamma_nu)*dG2 + dG3) + (- gamma_nu)^3*cfgaussder(time,nux,mu,sigma);

        case 4

            DG(iinu,:)  = (gamma_nu + 1)*((-gamma_nu)^3*dG1 + (-gamma_nu)^2*dG2 + (-gamma_nu)*dG3 + dG4) + (- gamma_nu)^4*cfgaussder(time,nux,mu,sigma);

        case 5

            DG(iinu,:)  = (gamma_nu + 1)*((-gamma_nu)^4*dG1 + (-gamma_nu)^3*dG2 + (-gamma_nu)^2*dG3 + (-gamma_nu)*dG4 + dG5) + (- gamma_nu)^5*cfgaussder(time,nux,mu,sigma);

    end

end

[T,N]   = meshgrid(time,nu);

%%

new_t       = time(dG1 <= 0);

ord_t       = linspace(new_t(1),time(end),numel(time));

dR          = df1dt1(new_t);

%% -- plot things

%

% parameters

graph.MarkerSize    = 10;

graph.LineWidth     = 1.5;

graph.MinorTicks    = 'on';

graph.FontName      = 'Serif';

graph.FontSize      = 30;

% plotting

imagesize            = [0.9 0.9];

papersize           = [60 20]; % in centimeters

figure('Color','w');%,'Units','Normalized','Position',[0.01 0.05 imagesize]);

set(gcf, 'PaperUnits','centimeters');

set(gcf, 'Units','centimeters');

set(gcf, 'PaperSize', papersize);

set(gcf, 'Position',[papersize.*(1 - imagesize)/2 imagesize.*papersize]);

set(gcf, 'PaperPositionMode', 'manual');

set(gcf, 'PaperPosition',[papersize.*(1 - imagesize)/2 imagesize.*papersize]);

set(gcf, 'Renderer', 'painters');

colours     = [jet(nu_nod);[0 0 0]];

dG_trad     = eval(sprintf('dG%d',n));

%Data = [DG./repmat(max(abs(DG),[],2),1,time_nod);dG_trad./max(abs(dG_trad))]; % normalised

Data  = [DGn;-sqrt(2*pi)*sigma*dR];

stars = [Par.Start,mu];

new_time = [new_time;new_t];

for ii = 1 : nu_nod+1,

    plot(new_time(ii,:)- stars(ii),Data(ii,:),...

        'LineWidth',graph.LineWidth,'Color',colours(ii,:)), %shading interp,

    hold on;

end

box on,

set(gca,'LineWidth',graph.LineWidth,'FontSize',graph.FontSize,'FontName',...

    graph.FontName,'XMinorTick',graph.MinorTicks,'YMinorTick',graph.MinorTicks,...

    'Units','Normalized','Position',[.07 0.1 0.88 0.84],...

    'TickLabelInterpreter', 'latex','YLim',[0 1]),

xlabel('$$\tau = t - \mu_{\nu}$$ (dimensionless)','FontSize',graph.FontSize,'FontName',graph.FontName,...

    'Interpreter','latex');

%ylabel('\nu-order','FontSize',graph.FontSize,'FontName',graph.FontName);

ylabel('$$R_{\sigma}^{\nu}(\tau)$$','FontSize',graph.FontSize,...

    'FontName',graph.FontName,'Interpreter','latex');

lgnd_   = strsplit([sprintf('$$\\nu = %.2f$$,',nu),...

    ['$$\nu = ',sprintf('%.2f',n),'$$,']],',');

lgnd    = lgnd_(1:end-1); xlim([0 5]),

legend(lgnd,'FontSize',graph.FontSize-10,'FontName',graph.FontName,...

    'Interpreter','Latex','Location','Best'),

childrengca = get(gca,'Children');

childrengca(1).LineStyle = '--';

tickValues = get(gca,'YTick');

newLabels = arrayfun(@(value)(sprintf('%.1f',value)), tickValues, 'UniformOutput',false);

set(gca, 'YTickLabel', newLabels);

tickValues = get(gca,'XTick');

newLabels = arrayfun(@(value)(sprintf('%.1f',value)), tickValues, 'UniformOutput',false);

set(gca, 'XTickLabel', newLabels);

🔗 参考文献

[1]卢少平、应黎明、王霞、李文杰.基于用户出行模拟的电动汽车快充站负荷预测及其优化调度[J].电力建设, 2020, 41(11):11.DOI:10.12204/j.issn.1000-7229.2020.11.004.

🍅往期回顾扫扫下方二维码