💥1 概述
在人类、工业和电动汽车的能源需求的推动下，全球能源需求预计将在未来几年稳步增长；更准确地说，预计到 2030 年增长将达到 40%。这种需求是由人类日益依赖能源的生活方式、电动汽车作为主要交通工具的出现以及机器将促进流程的进一步自动化。在当今的电网中，能源是在集中式和大型能源工厂（微电网能源发电）中生产的；然后，将能量传输到最终客户端，通常是通过非常远的距离并通过复杂的能量传输网格。如此复杂的结构，灵活性降低，难以适应需求增长，从而增加了电网不稳定和停电的概率。影响是巨大的，最近欧洲和北美的停电造成了数百万欧元的损失就证明了这一点。
在这种情况下，微电网正在成为一种有前途的能源解决方案，其中分布式（可再生）能源正在满足当地需求。当当地生产无法满足微电网要求时，从主要公用事业公司购买能源。设想微电网将提供许多好处，例如电力输送的可靠性（例如，通过岛屿），通过增加可再生能源的普及率来提高效率和可持续性，可扩展性和投资延期，以及提供辅助服务。孤岛是微电网的突出特征之一，是指能够断开微电网负荷与主电网的连接，并通过当地能源专门为其供电。在主电网无法支持总需求和/或运营商检测到一些可能退化为停电的主要电网问题的情况下，将执行预期孤岛。在这些情况下，微电网可以提供足够的能量，至少可以保证基本的电力服务。一旦整个系统再次稳定，将恢复与主电网的连接。显然，这些都是可能导致不稳定的重要功能。
为了提高智能电网的能力，一种典型的方法是考虑几个微电网相互交换能量的情况，即使这些微电网是孤岛的，即与主电网断开连接。换言之，在一组连续的微电网内部存在能量流，但在微电网和主电网之间不存在能量流动。在这种背景下，最优潮流问题最近引起了相当大的关注。例如，Ochoa 和 Harrison 联合考虑了协调电压控制的潮流问题。或者，相关文献中的工作侧重于不平衡配电网络，并提出了一种基于牛顿下降法解决三相潮流问题的方法。由于这些集中式解决方案可能会受到可扩展性和隐私问题的影响，一般来说，最优潮流问题是非凸的；因此，精确的解决方案可能太复杂而无法计算。因此，经常采用次优方法。一些文献中采用所谓的乘法器交替方向方法（参以分布式方式解决潮流问题。、
在本文中，为孤岛微电网之间的能源交易开发了一个分布式凸优化框架。更具体地说，该问题由几个孤岛微电网组成，这些电网通过任意拓扑交换能量流。由于可扩展性问题，为了保护成本函数的局部信息，讲解了一种基于次梯度的成本最小化算法，该算法在实际迭代次数内收敛到最优解。此外，这种方法允许进行非常直观的经济学解释，从“供需模型”和“市场清算”的角度解释算法迭代数值结果给出了算法的收敛速度和不同网络拓扑的实现成本。
📚2 运行结果
考虑一个由 M 个相互连接的微电网
组成的多微电网系统，该微电网在孤岛模式下运行。在每个调度间隔期间，每个微电网 
产生
个单位的电能并消耗
个单位的电能。此外，可以允许
向
出售电能
（
），并从
（
）购买能量
。然后，
内的电能平衡需要：
部分代码：
%% 创建成本函数 %=====发电机成本函数====== C = cell(1,M); C_prime = cell(1,M); C_primeInv = cell(1,M); for m=1:M [C{m}, C_prime{m}, C_primeInv{m}] = create_cost_func(generators{m}, .98*Pmax(m)); end %====变压器成本函数========== [gamma_cost, gamma_prime, gamma_primeInv] = create_cost_func(a, b, c, CableCapacity); %% 调试 if IsDebug || verbose lambda_fig = figure; title('Lambdas'); xlabel('迭代次数'); cost_fig = figure; title('Local Costs'); xlabel('iterations'); gap_fig = figure; title('Duality Gap'); xlabel('迭代次数'); ellipse = figure; hold on end %% 为 Lambda 找到（松散的）边界 Lambdas_min = inf; Lambdas_max = -inf; for mm = 1:M tmp = C_prime{mm}(0); if Lambdas_min > tmp Lambdas_min = tmp; end tmp = C_prime{mm}(1.01*max(Pmax(mm),E_c(mm))); if Lambdas_max < tmp Lambdas_max = tmp; end end % Initialize Lambdas and containing ellipsoid Lambdas = (Lambdas_min + Lambdas_max)/2*ones(M,1); radius2 = (Lambdas_max - Lambdas_min)^2/2; ElpsMatrix = eye(M) * radius2; if IsDebug || verbose updateplots(lambda_fig,1,Lambdas); if M == 2 figure(ellipse) line([Lambdas_min*[1 1] Lambdas_max*[1 1]], [Lambdas_min Lambdas_max*[1 1] Lambdas_min]) plotellipse(ElpsMatrix, Lambdas(1), Lambdas(2)) plot(Lambdas(1), Lambdas(2), 'xk') set(gca,'DataAspectRatio',[1 1 1]) end end %% 主要算法 E_min = zeros(M); E_sell = zeros(M,1); counter = 1; cost = zeros(M,1); cost2 = zeros(M,1); duality_gap = 10; Cbest = -inf; subg_nrmlz = inf; while counter < MaxIterations && (subg_nrmlz > MaxEpsilon)% || abs(duality_gap) > MaxDualityGap) counter = counter+1; %===对于给定的Lambdas，解决微电网局部问题==== for ii=1:M [E_sell(ii), E_min(:,ii)] = local_problem(E_c(ii), A(:,ii), Lambdas(ii), Lambdas, C_prime{ii}, C_primeInv{ii}, gamma_prime, gamma_primeInv); end %=====计算对偶成本函数的次梯度===== subgradient = sum(E_min,2) - E_sell; %====并对解椭球进行归一化处理====== subg_nrmlz = sqrt(subgradient'*ElpsMatrix*subgradient); subgradient = subgradient / subg_nrmlz; %====计算给定Lambdas的总成本===== tmp = 0; for mm = 1:M tmp = tmp + C{mm}(E_c(mm) + E_sell(mm) - sum(E_min(:,mm))) + sum(gamma_cost(E_min(:,mm))) + sum(Lambdas.*E_min(:,mm)) - Lambdas(mm)*E_sell(mm); end if Cbest < tmp Cbest = tmp; end %====更新Lambdas并包含椭球体============= alpha_step = 0;%(Cbest - tmp)/subg_nrmlz; jump = ElpsMatrix * subgradient; Lambdas_old = Lambdas; Lambdas = Lambdas + (1 + M*alpha_step)/(M+1) * jump; ElpsMatrix_old = ElpsMatrix; ElpsMatrix = M^2/(M^2-1)*(1-alpha_step^2)*(ElpsMatrix - 2*(1+M*alpha_step)/(M+1)/(1+alpha_step)*(jump*jump')); %====确保解决方案是可以接受的==== radi = eig(ElpsMatrix); if any(radi<-.1) alpha_step = 0; Lambdas = Lambdas_old + (1 + M*alpha_step)/(M+1) * jump; ElpsMatrix = M^2/(M^2-1)*(1-alpha_step^2)*(ElpsMatrix_old - 2*(1+M*alpha_step)/(M+1)/(1+alpha_step)*(jump*jump')); end radi = eig(ElpsMatrix); if any(radi)<-.1 fprintf('cucu') end %====检查新的 Lambas 是否在原来的范围内====== while any(Lambdas < Lambdas_min) for mm = 1:M if Lambdas(mm) < Lambdas_min subgradient = zeros(M,1); subgradient(mm) = -1; subg_nrmlz = sqrt(ElpsMatrix(mm,mm)); subgradient = subgradient / subg_nrmlz; alpha_step = (-Lambdas(mm)+Lambdas_min)/subg_nrmlz; jump = ElpsMatrix * subgradient; Lambdas = Lambdas - (1 + M*alpha_step)/(M+1) * jump; ElpsMatrix = M^2/(M^2-1)*(1-alpha_step^2)*(ElpsMatrix - 2*(1+M*alpha_step)/(M+1)/(1+alpha_step)*(jump*jump')); subg_nrmlz = inf; end end end %====这不完全是对偶差距，但有助于说明收敛性==== duality_gap = sum(Lambdas.*(sum(E_min,2)-E_sell)); if (IsDebug && mod(counter,1)==0) || (verbose && mod(counter,50)==0) fprintf('Iteration %d: volume = %e --- epsilon = %f --- duality gap = %f\n', counter, prod(radi), subg_nrmlz, duality_gap) fprintf('\nLambdas =\n'); fprintf(' %-10.4f\n',Lambdas); fprintf('E_min =\n'); fprintf([' ' repmat('%-10g',1,M) '\n'],E_min'); updateplots(lambda_fig,counter,Lambdas); fprintf('\nE_sell = \n'); fprintf(' %-10g \n',E_sell'); E_gen = E_c+sum(E_min,2)-sum(E_min,1)'; fprintf('\nE_gen = E_c = \n'); fprintf(' %-10g %-10g\n',[E_gen E_c]'); for m=1:M cost(m) = C{m}(E_gen(m)) + sum(gamma_cost(E_min(:,m))); cost2(m) = C{m}(E_c(m) + E_sell(m) - sum(E_min(:,m),1)) + sum(gamma_cost(E_min(:,m))) ... + sum(Lambdas.*E_min(:,m)) - Lambdas(m) * E_sell(m); end fprintf('\nTotal energy cost: %f (%f) USD per hour.\n\n', sum(cost), sum(cost2)); updateplots(cost_fig,counter,cost); updateplots(gap_fig,counter,duality_gap); if M==2 figure(ellipse) plotellipse(ElpsMatrix,Lambdas(1),Lambdas(2)) plot(Lambdas(1),Lambdas(2),'xk') end fprintf([repmat('-',1,50) '\n']); pause(1); end end if IsDebug updateplots(lambda_fig,counter,Lambdas); updateplots(cost_fig,counter,cost); updateplots(gap_fig,counter,duality_gap); end %% 输出 fprintf('\nSimulation ended after %d iterations.', counter); fprintf('\nFinal duality gap: %f.\n\n', duality_gap); fprintf('Lambdas =\n'); fprintf(' %-10.4f\n',Lambdas); fprintf('E_min =\n'); fprintf([' ' repmat('%-12g',1,M) '\n'],E_min'); fprintf('\nE_sell = \n'); fprintf(' %10g (%10g)\n',[ sum(E_min,2)'; E_sell']); E_gen = E_c+sum(E_min,2)-sum(E_min,1)'; fprintf('\nE_gen = E_c = \n'); fprintf(' %-12g %-12g\n',[E_gen E_c]'); gen_cost = zeros(M,1); tx_cost = zeros(M,1); for m=1:M gen_cost(m) = C{m}(E_gen(m)); tx_cost(m) = sum(gamma_cost(E_min(:,m))); cost2(m) = C{m}(E_c(m) + E_sell(m) - sum(E_min(:,m),1)) + sum(gamma_cost(E_min(:,m))) ... + sum(Lambdas.*E_min(:,m)) - Lambdas(m) * E_sell(m); end cost = gen_cost + tx_cost; uGcosts(:,icase) = cost2; fprintf('\nTotal energy cost per microgrid:') fprintf('\n Generation Transmission Total\n') fprintf(' %-12g %-12g %-12g\n',[gen_cost tx_cost cost]') fprintf('\nTotal energy cost: %f (%f) USD per hour.\n\n', sum(cost), sum(cost2)); Sell4(icase) = sum(E_min(M,:)); Price(icase) = Lambdas(M); end %% 可视化 figure plot(inputdata{:,end-2},uGcosts'/1000) xlabel('E_4^{(c)} [MWh]') ylabel('成本[k$]') grid axis([0 12 0 1.4]) hold plot(inputdata{:,end-2},C{1}(inputdata{:,1})/1000,'k--') plot(inputdata{:,end-2},C{1}(inputdata{:,end-2})/1000,'k-.') legend('uG1','uG2','uG3','uG4','disc. {1,2,3}','disc. 4','Location','northwest') saveas(gcf, 'local_costs.png') figure hh = plotyy(inputdata{:,end-2},Sell4,inputdata{:,end-2},[Price Sell4.*Price]/1000); xlabel('E_4^{(c)} [MWh]') ylabel(hh(1),'[MWh]') ylabel(hh(2),'[k$]') grid axis(hh(1),[0 12 0 3.6]) set(hh(1),'YTick',0:3) axis(hh(2),[0 12 0 1.8]) set(hh(2),'YTick',0:.5:1.5) legend('售出电能', '单价', '收益') saveas(gcf, 'trading.png')
🎉3 参考文献
部分理论来源于网络，如有侵权请联系删除。
[1]David Gregoratti, Javier Matamoros (2018) Distributed Energy Trading: The Multiple-Microgrid Case.
[2]李长宇,唐文秀.基于数据驱动的多微电网互联系统分布鲁棒运行优化[J].智慧电力,2022,50(05):1-8.
