MATLAB代码:考虑电动汽车时空调度的配电网潮流优化 参考文档:He L , Yang J , Yan J , et al. A bi-layer optimization based temporal and spatial scheduling for large-scale electric vehicles[J]. Applied Energy, 2016, 168(apr.15):179-192. DOI:10.1016/j.apenergy.2016.01.089 主要内容:机组组合采用原文相同的线性化方法最优潮流采用二阶锥松弛替代原文算法前阵子帮小区物业调配电系统物业大叔愁得头发都掉了几根——夏天晚上七八点大家下班回家开空调、充电动车配电箱直接过载跳闸连路灯都暗了半截。其实这就是电动汽车并网调度最常见的痛点用户想省钱找谷时段充电但扎堆充电又会把电网搞崩得有个法子平衡两边的需求。MATLAB代码:考虑电动汽车时空调度的配电网潮流优化 参考文档:He L , Yang J , Yan J , et al. A bi-layer optimization based temporal and spatial scheduling for large-scale electric vehicles[J]. Applied Energy, 2016, 168(apr.15):179-192. DOI:10.1016/j.apenergy.2016.01.089 主要内容:机组组合采用原文相同的线性化方法最优潮流采用二阶锥松弛替代原文算法刚好翻到2016年那篇发在Applied Energy上的经典论文双层调度的思路挺绝但原文用的下层迭代算法算起来太磨人我这次换了二阶锥松弛来做最优潮流再配上和原文一致的机组组合线性化套路还写了个能直接跑的Matlab代码今天就把踩过的坑和干货都唠出来。先唠明白为啥要这么改原文的双层框架其实很好懂上层管配电网的机组启停和出力计划下层管每台EV车主的充电需求两边互相妥协找最优解。但原文下层用的非线性潮流求解不仅慢还容易卡在局部最优换成二阶锥松弛之后就能把非线性的潮流方程转换成凸优化问题求解器几秒就能出全局最优解而且机组组合的线性化套路和原文完全对齐不用改太多逻辑。直接上代码边写边拆咱们用最简单的3节点配电网当测试案例节点1是主网联入的Slack节点节点2和3各装了5台充电桩调度周期是24小时每小时一个时段。首先初始化参数和工具箱%% 考虑EV调度的配电网潮流优化二阶锥松弛版 % 对齐Applied Energy 2016论文的机组组合逻辑替换下层潮流为二阶锥松弛 clear;clc;close all % 装了YALMIP和Gurobi的同学直接跑没装的换intlinprog也能跑就是慢第一步初始化系统参数先把配电网、发电机、EV的基础参数都写清楚不用太复杂能跑通就行%% 1. 系统参数初始化 bus_num 3; % 总节点数 gen_num 1; % 配网内传统发电机组数量 ev_num_per_bus 5; % 每个充电桩节点的EV数量 T 24; % 调度时段24小时 % 发电机组参数 P_gen_max 10; % 最大出力MW P_gen_min 2; % 最小出力MW R_up 3; R_down 3; % 上下爬坡速率MW/h fuel_cost 0.1; % 发电燃料成本 元/MWh % EV参数 ev_cap 50; % 单台EV电池容量kWh ev_charge_max 3; % 单台EV最大充电功率kW % 车主默认下班回家要把电池从20%充到90%单台总充电需求35kWh ev_total_demand 35*ev_num_per_bus/1000; % 转成MW·h % 基础负荷不含EV充电 base_load [5;3;4]; % 节点1是Slack节点2、3的基础负荷 % 配网支路参数支路1-21-3阻抗rjx branch [0.01 0.1; 0.02 0.15]; bus_from [1;1]; bus_to [2;3]; % 峰谷电价白天8-22点电价高 price ones(1,T)*0.5; price(8:22) 1.2;第二步定义优化变量这里用YALMIP的变量语法比自己手写矩阵方便太多%% 2. 定义优化变量 % 发电机启停状态二进制变量和出力 u_gen binvar(gen_num,T); P_gen sdpvar(gen_num,T); % 二阶锥松弛专用变量节点电压平方避开极坐标的非线性sin/cos V2 sdpvar(bus_num,T); % EV充电功率节点2、3的总充电功率 P_ev_total sdpvar(2,T); % 行对应节点2、3列对应时段 % 支路有功/无功功率用来写潮流约束 P_branch sdpvar(size(branch,1),T); Q_branch sdpvar(size(branch,1),T);第三步添加约束和原文对齐首先是机组组合的线性化约束和原文的逻辑完全一致%% 3. 添加约束 Constraints []; % -------------------------- 机组组合约束 -------------------------- for t 1:T % 发电机出力必须在启停状态对应的上下限之间 Constraints [Constraints, P_gen(:,t) P_gen_min*u_gen(:,t)]; Constraints [Constraints, P_gen(:,t) P_gen_max*u_gen(:,t)]; % 爬坡约束防止机组出力突变 if t1 Constraints [Constraints, P_gen(:,t)-P_gen(:,t-1) R_up*u_gen(:,t) P_gen_max*(1-u_gen(:,t-1))]; Constraints [Constraints, P_gen(:,t-1)-P_gen(:,t) R_down*u_gen(:,t-1) P_gen_max*(1-u_gen(:,t))]; end end % 初始状态默认机组开机避免凌晨突然停机 Constraints [Constraints, u_gen(:,1)1]; % -------------------------- EV充电约束 -------------------------- for b1:2 % 只处理节点2、3的EV for t1:T % 单节点总充电功率不能超过上限 Constraints [Constraints, P_ev_total(b,t) ev_num_per_bus*ev_charge_max/1000]; Constraints [Constraints, P_ev_total(b,t) 0]; % 全天总充电量必须满足车主需求 if tT Constraints [Constraints, sum(P_ev_total(b,:)) ev_total_demand]; end end end % -------------------------- 配网潮流约束二阶锥松弛版 -------------------------- % Slack节点电压固定为1.0p.u.平方就是1 Constraints [Constraints, V2(1,:)1]; % 节点电压必须在合格范围内0.95~1.05p.u. for t1:T Constraints [Constraints, 0.95^2 V2(:,t) 1.05^2]; end % 支路功率的二阶锥约束还有节点有功平衡 for k1:size(branch,1) % 计算支路导纳 r branch(k,1);xbranch(k,2); g r/(r^2x^2); b -x/(r^2x^2); for t1:T Vi sqrt(V2(bus_from(k),t)); Vj sqrt(V2(bus_to(k),t)); % 二阶锥松弛核心支路功率模长不超过电压平方的组合 Constraints [Constraints, norm([P_branch(k,t);Q_branch(k,t)]) sqrt(V2(bus_from(k),t)*V2(bus_to(k),t))*sqrt(g^2 b^2)]; % 节点有功平衡总注入功率 基础负荷 EV负荷 - 支路流出功率 for b_node2:bus_num if b_node 1 net_power P_gen(1,t) - base_load(b_node) - P_ev_total(b_node-1,t); else net_power -base_load(b_node) - P_ev_total(b_node-1,t); end % 统计所有连接到该节点的支路功率 branch_flow 0; if bus_from(k)b_node branch_flow branch_flow P_branch(k,t); end if bus_to(k)b_node branch_flow branch_flow - P_branch(k,t); end Constraints [Constraints, net_power branch_flow]; end end end第四步定义目标函数并求解目标函数和原文一致最小化总运行成本包括发电燃料成本和EV充电的峰谷电价成本%% 4. 目标函数与求解 % 总发电成本 gen_cost fuel_cost * sum(sum(P_gen)); % 总EV充电成本MW·h * 元/MWh 元 ev_cost sum(sum(P_ev_total.*repmat(price,2,1))); Objective gen_cost ev_cost; % 调用求解器用Gurobi最快没装的换solver,intlinprog ops sdpsettings(solver,gurobi,verbose,1); result optimize(Constraints,Objective,ops);第五步结果展示跑通之后直接画图看结果%% 5. 结果展示 if result.problem 0 disp(恭喜求解成功啦); P_gen_res value(P_gen); V2_res value(V2); P_ev_res value(P_ev_total); % 画发电机出力曲线 figure(Name,机组出力); plot(1:T,P_gen_res(1,:),linewidth,2); xlabel(时段(h));ylabel(出力(MW));title(传统机组出力);grid on; % 画EV充电曲线 figure(Name,EV充电调度); plot(1:T,P_ev_res(1,:),r,1:T,P_ev_res(2,:),b,linewidth,2); xlabel(时段(h));ylabel(总充电功率(MW));legend(节点2,节点3);grid on; % 画节点电压 figure(Name,节点电压); plot(1:T,sqrt(V2_res(2,:)),r,1:T,sqrt(V2_res(3,:)),b,1:T,ones(1,T),k--,linewidth,1.5); xlabel(时段(h));ylabel(电压(p.u.));legend(节点2,节点3,标准电压);grid on; else disp(坏啦求解失败错误码);disp(result.problem); end踩过的坑唠两句一开始我直接用非线性潮流写结果求解器跑了十分钟都没出结果换成二阶锥松弛之后秒出。不过要注意几个坑二阶锥松弛对于径向配电网是精确的不用担心里解不出来可行解但如果是网状配网可能需要额外的约束一开始忘了加电压上下限结果算出来电压跑到1.5p.u.明显不符合电网标准EV的充电约束一开始只加了单台功率上限忘了总充电量要满足车主的需求结果算出来EV根本没充满电。最后补两句这个代码只是个入门demo要是真的做大规模配网调度还得加上EV的到达/离开时间、储能装置、不同充电桩的功率限制这些细节但作为理解双层调度和二阶锥松弛的例子已经够用了。要是把代码改成和原文完全一致的双层迭代结构也只需要把下层的EV调度单独抽出来循环求解就行感兴趣的同学可以自己改改~