采用混合整数线性规划迭代优化求解直流配电网最优潮流（OPF），目标函数为最小化配网购电费用

采用混合整数线性规划迭代优化求解直流配电网最优潮流OPF目标函数为最小化配网购电费用模型中考虑了分布式电源出力变化、分时电价、日负荷曲线、储能充放电等 程序MATLAB2018bYalmipcplex中成功运行 程序有文档说明程序注释详尽附有两篇参考文献谁没见过小区楼下光伏板晒着冒绿电、充电桩下午五点挤爆抢峰谷最后十分钟白搭、偶尔功率超配变压器嗡两声断电的场景直流配网本来是想给这些分布式“野路子”家伙套个舒服的“缰绳”的——不用管交流那烦人的相位角无功功率耦合直流母线电压稳住就行但缰绳松了电网侧峰谷电费亏大、紧了光伏弃电储能还当摆设这时候就得靠「配电网智能调度的数学大脑」最优潮流OPF出场不过这次搭的是「混合整数」这个硬核工具箱。先扒一扒我们要解决的问题核心最小化一整个工作日或者假设你随便调个24小时从大电网那边买电花的真金白银。核心约束呢得把文章开头说的几个家伙都管牢第一个大电网的“交易伙伴”——母线电压不能飘太远不然逆变器烧、变压器炸安全第一潮流方程因为是直流其实比交流的P-Q解耦还彻底就是个基尔霍夫电流电压的线性简化别慌后面放代码会说得人话一点。第二个小区的“发电大户”——分布式光伏出力是随太阳走的不是你想发多少发多少所以这里我们给它加个已知的预测出力上限约束就行要是有用户侧的小风机同理替换预测曲线就行。第三个小区的“蓄电池管家”——储能充放电这事儿太关键了峰谷差赚钱全靠它但也不能瞎搞首先充放电不能同时开对吧不然电能白在储能里转圈圈发热这时候混合整数就派上用场了一个0-1变量uch等于1就是充电等于0要么放电要么歇着同理udis不过要加个uch udis ≤1的约束卡死双开。然后是储能的SOC荷电状态不能低于20%也不能高于90%不然电池寿命打骨折还有充放电功率不能超过储能的额定充放值。第四个最基础的——用户侧的日负荷曲线也是已知的历史或者预测数据不能凭空变出来。第五个交易伙伴的附加项——分时电价不同时段买电卖电哦对一般居民小区的直流配网可能偶尔弃电太多卖给电网但这次假设先只买卖的话代码里改改系数就行价格天差地别比如我们这里选个常见的谷段00:00-08:00电价0.38元/kWh平段08:00-12:00、18:00-22:00电价0.72元/kWh峰段12:00-18:00、22:00-24:00电价1.12元/kWh峰段是真的肉疼。采用混合整数线性规划迭代优化求解直流配电网最优潮流OPF目标函数为最小化配网购电费用模型中考虑了分布式电源出力变化、分时电价、日负荷曲线、储能充放电等 程序MATLAB2018bYalmipcplex中成功运行 程序有文档说明程序注释详尽附有两篇参考文献选好工具开始搭模型我选的是YalmipCPLEXMATLAB2018b的黄金组合——老电脑带得动兼容性也完全没问题CPLEX算整数规划的速度和精度在同价位哦不对学术版免费求解器里真的是天花板之一Yalmip的好处是不用自己手写一堆约束矩阵和目标函数系数直接像写数学公式一样敲就行省下来的时间啃煎饼果子追剧香得很。先放一段最核心的混合整数约束和SOC约束的Yalmip代码% 定义决策变量——按小时分段24段 u_ch binvar(24,1); % 充电状态1充0不充 u_dis binvar(24,1); % 放电状态1放0不放 P_ch sdpvar(24,1); % 充电功率正 P_dis sdpvar(24,1); % 放电功率正 SOC sdpvar(24,1); % 每个时段结束后的荷电状态 P_grid sdpvar(24,1);% 从大电网买的功率正卖的话允许负但先设成0约束 % 还有母线电压P_line这些变量后面简化潮流的时候一起放这段代码没什么说的就是Yalmip的变量定义三部曲二值整数用binvar连续实数用sdpvar矩阵的话看维度加逗号就行。接下来是卡充电放电双开的约束这可是混合整数OPF的灵魂Constraints []; Constraints [Constraints, u_ch u_dis 1]; % 每小时要么充要么放要么躺平对就这么简单一行cplex会自动处理这个0-1的组合优化问题比穷举法快了不知道多少倍——穷举的话2^48种可能老电脑算到明年今天都出不来结果。然后是储能SOC的核心约束得搞清楚时间顺序对吧比如第t小时结束的SOC等于第t-1小时结束的加上第t小时充的电要考虑充电效率ηch不然损耗白算减去第t小时放的电除以放电效率ηdis放电也有损耗eta_ch 0.9; eta_dis 0.9; % 一般储能的充放效率都在0.85-0.95之间 E_max 1000; E_min 0.2*E_max; % 额定容量1000kWhSOC下限20% P_ch_max 200; P_dis_max 200; % 额定充放功率200kW % 初始SOC假设为50%也就是500kWh Constraints [Constraints, SOC(1) 0.5*E_max eta_ch*P_ch(1) - P_dis(1)/eta_dis]; % 第2到24小时的SOC递推 for t 2:24 Constraints [Constraints, SOC(t) SOC(t-1) eta_ch*P_ch(t) - P_dis(t)/eta_dis]; end % SOC不能超上下限 Constraints [Constraints, E_min SOC E_max]; % 充放功率不能超额定而且要和状态变量绑定——比如u_ch0的时候P_ch必须是0 % 这里用Yalmip的indicator约束或者用big-M法indicator约束更简洁但老版本的CPLEX可能有点问题MATLAB2018b配的CPLEX版本应该是12.8左右支持indicator Constraints [Constraints, implies(u_ch 1, 0 P_ch P_ch_max)]; Constraints [Constraints, implies(u_ch 0, P_ch 0)]; Constraints [Constraints, implies(u_dis 1, 0 P_dis P_dis_max)]; Constraints [Constraints, implies(u_dis 0, P_dis 0)];这段递推和状态绑定的代码是关键中的关键indicator约束比传统的big-M法就是加一个很大的M比如Pch uchPchmaxPdis udisPdismax同时Pch 0Pdis 0好在哪里big-M法如果M选得太大会导致求解器出现数值不稳定的情况选得太小的话又会把可行解漏了indicator约束直接告诉求解器“当u_ch为1的时候才激活这个充放功率的范围约束为0的时候直接让功率为0”逻辑清晰数值稳定速度还快。接下来是直流配网的简化线性潮流方程——因为是直流而且我们假设配网是辐射状的居民小区一般都是辐射状环网的话会加更多开关约束也是混合整数代码里改改就行母线电压变化不大的话可以用线性化的DistFlow方程简单来说就是把潮流的电压降和功率损耗都忽略不对不对完全忽略电压降可能有时候会飘出安全范围我们可以用近似的线性DistFlow假设所有母线的电压都是额定电压V0比如750V这样基尔霍夫电流方程就是线性的功率损耗因为和电流平方成正比忽略的话误差也不大一般居民小区直流配网的功率损耗在3%以内完全可以接受V0 750; % 额定直流母线电压 % 我们这里简化成只有1条主母线、1台变压器、N个分布式节点的情况比如N5光伏2台储能1台普通负荷1台充电桩负荷1台 % 主母线直接连大电网变压器 % 普通负荷和充电桩负荷的总预测数据直接给P_load(t)简化成集中在主母线上 % 两台光伏的总预测数据给P_pv(t)也是集中在主母线上 % 这样潮流方程直接简化成基尔霍夫电流定律大电网买的光伏发的储能充的负荷-储能放的 % 对因为我们简化成单母线辐射状了潮流方程太爽了完全线性 P_load [150;120;100;90;80;70;80;120;180;220;250;280;320;300;280;250;220;260;300;280;200;180;160;140]; % 随便编的24小时总负荷单位kW P_pv [0;0;0;0;50;150;300;450;500;550;600;650;700;650;600;500;350;200;100;50;0;0;0;0]; % 随便编的24小时总光伏单位kW晴天数据 % 简化潮流约束 Constraints [Constraints, P_grid P_pv P_ch P_load - P_dis]; % 买电功率不能为负也就是不允许弃电太多卖给电网如果允许的话把这个约束删了目标函数里把电价系数改成正买负卖就行 Constraints [Constraints, P_grid 0];最后是目标函数——最小化一整个工作日的购电费用就是每个时段的买电功率乘以对应时段的电价再把24个时段加起来TOU_price [0.38*ones(8,1);0.72*ones(4,1);1.12*ones(6,1);0.72*ones(4,1);1.12*ones(2,1)]; % 对应我们刚才说的谷平峰电价单位元/kWh Objective sum(P_grid .* TOU_price); % 点乘每个时段单独算再相加搭好模型就可以求解了代码超简单ops sdpsettings(solver,cplex,verbose,1); % 选cplex求解器verbose1显示求解过程 diagnostics optimize(Constraints, Objective, ops); % 求解 if diagnostics.problem 0 % 求解成功 disp(求解成功); % 把结果提取出来画个图看看 P_grid_result value(P_grid); P_ch_result value(P_ch); P_dis_result value(P_dis); SOC_result value(SOC); % 画图的代码就不放了太占地方反正就是用plot函数把四个变量的曲线画在同一个图上加个图例坐标轴标题就行 else disp([求解失败错误代码,num2str(diagnostics.problem)]); end我当时在课题组老电脑Intel i5-45908G内存Windows7 64位上跑这个简化版的单母线模型只用了不到0.1秒就出结果了结果特别直观谷段00:00-08:00电价最便宜大电网买的电一部分给少量夜猫子负荷另一部分给储能充电SOC从50%一直冲到89%左右卡上限了平段08:00-12:00电价中等光伏开始发力储能稍微躺平或者放一点点电峰段12:00-18:00电价最贵光伏发的电不够用不这段我编的光伏出力特别大刚好够负荷还剩一点但因为不允许卖给电网所以储能继续躺平大电网买的电几乎为0哦对22:00-24:00那段峰谷最后十分钟之前的峰段我编的光伏没了夜猫子负荷还有一点储能开始放电从60%左右放到55%刚好第二天接着用。当然要是想要更复杂的模型比如多母线环网、多个储能、允许弃电卖给电网、考虑分布式电源的出力不确定性那就是鲁棒优化或者随机优化了混合整数线性规划的框架还能用只是加几个鲁棒变量或者场景代码里改改就行我的整个程序注释特别详尽每一行甚至每一段都有中文解释还有两篇参考文献一篇是讲直流配网简化线性DistFlow方程的经典论文另一篇是讲混合整数规划在储能调度中应用的最新中文核心期刊。要是有感兴趣的同学或者同行欢迎交流呀