1主要内容
该程序复现《面向配电网韧性提升的移动储能预布局与动态调度策略》，具体摘要内容见下图，程序主要分为两大模块，第一部分是灾前预防代码，该部分采用两阶段优化算法，以移动储能配置成本和负荷削减成本作为目标函数，考虑光伏出力的不确定性，采用大M和列和约束生成算法（CC&G）；第二部分是灾后恢复代码，该部分未采用两阶段鲁棒，采用混合二阶锥规划算法进行优化，可求解得到最优负荷削减和设备出力方案。该程序采用gurobi作为求解器，kkt条件采用matlab函数直接求解，程序注释清晰。（该程序在网上比较流行，但是还是存在几个问题，下面和大家一一讨论）
- 灾前预防阶段
灾前预防采用CCG算法，运行一下程序得到的结果如下：
为了更加清晰看到CC&G算法收敛情况以及不确定变量的取值情况，补充一下两个结果图。
通过上述两个结果图看出，程序是两次即完全收敛，第二个图是不确定性分布式光伏出力（仅考虑5个时刻），从图中能够看到光伏除第一个时间点受鲁棒保守度限制为下限外，其他时刻均为上限，而目标为储能配置和负荷削减成本，这样的结果和常理不符，光伏出力越小应该成本会越高，为了进一步验证程序问题，将光伏设置到下限值。
从预布局成本结果来看，该目标值更大了，结合如下目标能够看出，光伏作为不确定变量，应该极大提升整体成本，从而给出最鲁棒的结果，通过该结果给运营商作为参考，检验最恶劣情况下的成本是否在接受范围。
但是通过这次测试，发现优化方向反了，导致结果不对，程序超低价分享给大家进一步研究，可以改变不确定度再检验一下，该种情况的具体问题可以通过检验对偶性来分析，我会抽空再做一个介绍。
- 灾后恢复阶段
灾后恢复由于变量数量比较多，运行时间比较长，matlab程序里面内置了选项，可以选择直接通过历史数据来形成结果图（如果不选历史数据，采用运行的方式，程序代码必须 改动一下才能运行成功，根据电脑配置运行时间存在差异，基本在十个小时左右）。
2部分代码
% 两阶段鲁棒优化问题，考虑移动储能预配置 %% 清理内存空间 clc clear close all warning off % % 随机生成初始场景it_max=5;%最大迭代次数 P_PV = zeros(5,it_max+1); P_PV(:,1) = [120 240 360 360 180]'/10000; P_PV0 = [120 240 360 360 180]/1000; % 光伏预测出力，MW tau = 0.45; % 不确定度 mpc = IEEE33; % 配电网参数 SB = mpc.baseMVA; % 基准功率,MW %% 迭代求解两阶段鲁棒优化问题 LB=zeros(1,it_max); UB=zeros(1,it_max); e=zeros(1,it_max); for it=1:it_max yalmip(' clear ') [alpha_ME , alpha_ij , LB(it)]=Master_Problem(P_PV(:,1:it),it); [P_PV(:,it+1),UB(it)]=Sub_Problem(alpha_ME , alpha_ij); if it>1 UB(it)=min(UB(it-1:it)); end e(it)=abs(UB(it)-LB(it))/UB(it); if abs(e(it))<1e-5 LB(it+1:it_max) = LB(it); UB(it+1:it_max) = UB(it); break end end %% 输出结果 disp(' ***************************两阶段鲁棒优化结果(考虑移动储能预配置)*************************** ') disp([' 断开的支路： ',num2str(find(~alpha_ij)' )])disp(['移动储能预配置节点：',num2str(find(alpha_ME)')]) Cpre = value(UB(it)); disp([' 预布局成本： ' , num2str(Cpre) , ' 元 ']) mpc = IEEE33; % 配电网参数 O_L0 = mpc.branch(:,1:2); % 支路集合 O_L0(alpha_ij == 0,:)=[]; O_DG = [2 9 29 6 25 13 18 20 24 33]; % DG节点集合，1-5为DEG,6-10为PV O_PV = O_DG(6:10); % PV节点集合 figure %节点的位置 node_location=[(1:14)' , 2 *ones( 14 , 1 );( 15 : 18 ) ',2.5*ones(4,1);(2:5)' ,ones( 4 , 1 );( 3 : 5 ) ',4*ones(3,1);(6:13)' , 3 *ones( 8 , 1 )];for k=1:length(O_L0) m = [node_location(O_L0(k,1),1),node_location(O_L0(k,2),1)]; n = [node_location(O_L0(k,1),2),node_location(O_L0(k,2),2)]; plot(m , n , 'k-' , 'LineWidth' , 1) hold on end axis off title('两阶段鲁棒优化的移动储能预配置与配电网重构结果'); essential_user = [3 4 6 10 11 15 17 19 24 26 28 33]; ordinary_user = setdiff(1:33,essential_user); h0 = scatter(node_location(essential_user,1) , node_location(essential_user,2) , 'ko' ,'filled' , 'LineWidth' , 1 ); h1 = scatter(node_location(ordinary_user,1) , node_location(ordinary_user,2) , 'ko' , 'LineWidth' , 1); h2 = scatter(node_location(O_DG(1:5),1) , node_location(O_DG(1:5),2) , 'r*' , 'SizeData' , 50 , 'LineWidth' , 1); h3 = scatter(node_location(O_PV,1) , node_location(O_PV,2) , 'rsquare' , 'SizeData' , 50 , 'LineWidth' , 1); h4 = scatter(node_location(alpha_ME == 1,1) , node_location(alpha_ME == 1,2) , 'r^' , 'SizeData' , 50 , 'LineWidth' , 1); legend([h0 h1 h2 h3 h4],{'重要用户' , '非重要用户' , 'DEG节点' , 'PV节点' , '移动储能预配置节点'}) figure; plot(UB-LB,'b-','LineWidth',1.5) xlabel('迭代次数'); ylabel('UB-LB'); figure; plot(P_PV0'/SB,'b-','LineWidth',1.5); hold on plot(P_PV0' /SB-P_PV 0 '/SB*tau,' k- ',' LineWidth ',1.5); plot(P_PV(:,3),' r-- ',' LineWidth ',1.5); legend(' 上限 ',' 下限 ',' 结果 '); xlabel(' 时间 '); ylabel(' 功率 ');
