yalmip+cplex机组组合调度代码实现案例

39-机组组合调度代码yalmipcplex 摘要代码是一个用来学习机组组合、熟悉yalmipcplex编程的好案例按照机组组合的约束以及难度依次设置了四个不同的case案例1为具有线性开关的简单机组组合案例2为在案例1的基础上添加了最小开关机时间约束案例3在2的基础上加入了量化功率平衡案例四则在案例3基础上实现了多天多次循环模拟 4、实现效果见下图搞电力系统入门仿真的谁没踩过纯写混合整数线性规划的坑啊变量加个下标i就忘到下个约束外开关机时间逻辑卡半天算不出可行解解出来又不知道是不是最优的。还好摸到了YALMIP搭模型配CPLEX算的路子不用自己凑约束矩阵爽到飞起。刚好手边有个分四个难度的机组组合入门练手代码拆解开唠唠新手照着敲绝对能通老手也能捡捡细节。先给代码搭个基础架子吧——先定义全局或者说通用的参数省得每个case都改。比如整个模拟的小时数、天数case4要用、机组的台数还有每台机组的参数最大最小出力、启停成本、线性的发电成本、最小开机停机时间这些。参数矩阵用cell或者struct都可以这里为了新手看着顺用二维数组存机组参数一维数组存负荷。% 通用参数 nt 24; % 单个周期24小时 case4_days 3; % case4的循环天数 ng 3; % 3台机组够演示了多了改参数就行 load [200 220 240 260 280 300 350 400 380 360 340 320 310 300 320 340 360 380 400 380 350 300 280 250]; % 典型日负荷曲线 % 机组参数矩阵行是机组i列是Pmin Pmax C_start C_stop a(线性成本系数) T_min_up T_min_down gen_params [ 50 300 1000 500 0.5 4 2; 100 400 2000 1000 0.4 6 4; 200 500 4000 2000 0.3 8 6; ];接下来先啃最简单的case1带线性开关约束的机组组合。啥是线性开关约束就是别让机组同时既开机又停机还有出力在停机时必须是0开机时在Pmin和Pmax之间——这些都能用YALMIP的整数变量和连续变量搭起来。先声明变量连续变量P(i,t)第i台机t时刻的出力二进制变量u(i,t)第i台机t时刻的状态1开0关哦对还有启停变量但case1里虽然可以算成本不过约束里可以先不用太纠结提前时间的问题不对等下开关约束本身还是要防止跳变太奇怪吗不case1标题说的是“具有线性开关”应该是指最最基础的出力上下限由u控制启停状态也可以表示成u的前后差但暂时不用加最小时间的硬约束对吧对不然就是case2了。% --- Case1: 带线性开关的简单UC --- yalmip(clear); % 每次case前清一下YALMIP的变量避免串 P sdpvar(ng, nt, full); % 连续出力 u binvar(ng, nt, full); % 状态二进制变量 cost 0; % 初始化总成本 % 目标函数总发电成本启停成本 for t 1:nt for i 1:ng cost cost gen_params(i,5)*P(i,t); % 线性发电成本 if t 1 % 开机成本t-1关t开即u(i,t)-u(i,t-1)1 cost cost gen_params(i,3)*max(0, u(i,t)-u(i,t-1)); % 关机成本t-1开t关同理 cost cost gen_params(i,4)*max(0, u(i,t-1)-u(i,t)); end end end % 约束集合 Constraints []; % 1. 功率平衡约束所有开着的机的出力加起来等于负荷case1先不考虑旋转备用简单点 for t 1:nt Constraints [Constraints, sum(P(:,t)) load(t)]; end % 2. 线性开关约束核心 for t 1:nt for i 1:ng % 停机时出力必须是0开机时在Pmin-Pmax之间 Constraints [Constraints, P(i,t) gen_params(i,1)*u(i,t)]; Constraints [Constraints, P(i,t) gen_params(i,2)*u(i,t)]; end end % 求解器设置选cplex算快准狠的混合整数规划 ops sdpsettings(solver,cplex,verbose,1); % verbose1看求解过程新手可以开着后面嫌烦改成0 result optimize(Constraints, cost, ops); % 简单看一下结果比如总成本还有u矩阵和P矩阵 if result.problem 0 disp([Case1 最优总成本, num2str(value(cost)), 元]); disp(Case1 机组状态矩阵行机组列小时1开0关); disp(value(u)); else disp(Case1 求解失败); end这里插一句线性开关约束的小细节——别小看这两个不等式新手一开始容易写成“如果u1则P在中间否则P0”这种逻辑判断YALMIP虽然能做逻辑约束比如用implies函数但implies有时候会转化成更复杂的约束或者引入额外变量这种用线性乘法的形式虽然u是二进制本质上是分段线性的整数约束写法更高效CPLEX最喜欢这种规规矩矩的MILP约束了。比如u0的时候两个不等式直接把P夹死在0u1的时候就变成普通的出力上下限完美接下来升级到case2加最小开关机时间约束。这个是电力系统机组组合里最常碰到的硬约束之一——火电机组不能说开就开说停就停不然锅炉、汽轮机寿命会折损得很快。怎么用YALMIP写这个约束呢首先我们要回忆一下最小开关机时间的逻辑最小开机时间Tminup如果一台机组在t时刻开机了也就是u(i,t)1且u(i,t-1)0那么接下来的Tminup-1个小时也就是t1到tTminup-1它必须保持开机状态不能停。最小停机时间Tmindown同理如果t时刻停机了u(i,t)0且u(i,t-1)1接下来Tmindown-1个小时必须保持停机。但逻辑判断直接写implies的话刚才说了可能不够爽或者试试更直观的前后差变量我们可以再声明两个二进制变量s_up(i,t)开机动作1表示t时刻开机0否则等于max(0, u(i,t)-u(i,t-1))不过用binvar声明的话约束起来更直接s_down(i,t)同理停机动作不过其实用u本身加上循环的边界条件处理t1的时候假设是什么状态这里新手入门可以假设所有机组在模拟开始前的Tminup-1和Tmindown-1小时内都是停机状态——这样第一台开机的机组不用考虑前序时间的限制不过实际工程里可以设更合理的初始状态也能写出来。% --- Case2: Case1基础上加最小开关机时间 --- yalmip(clear); P sdpvar(ng, nt, full); u binvar(ng, nt, full); cost 0; % 目标函数和case1一样 for t 1:nt for i 1:ng cost cost gen_params(i,5)*P(i,t); if t 1 cost cost gen_params(i,3)*max(0, u(i,t)-u(i,t-1)); cost cost gen_params(i,4)*max(0, u(i,t-1)-u(i,t)); end end end Constraints []; % 功率平衡和线性开关约束照搬case1 for t 1:nt Constraints [Constraints, sum(P(:,t)) load(t)]; for i 1:ng Constraints [Constraints, P(i,t) gen_params(i,1)*u(i,t)]; Constraints [Constraints, P(i,t) gen_params(i,2)*u(i,t)]; end end % --- 新增最小开关机时间约束 --- % 先处理初始边界条件假设t1时所有机组都停机 for i 1:ng % 最小开机时间s_up(i,t)触发后接下来T_min_up-1小时u1 % 从t1开始算先找开机动作可能影响的时间段 for t 1:nt if t gen_params(i,6) - 1 nt continue; % 最后几个小时开机的话不需要强制到模拟结束外 end % 用implies其实新手更易懂直接翻译逻辑如果t时刻开机那么t到tT-1都是1 Constraints [Constraints, implies(u(i,t)-u(i,t-1)1, u(i,t:tgen_params(i,6)-1)1)]; % 这里要注意t1的时候u(i,0)不存在哦刚才漏写了t1的边界implies % 补t1的最小开机约束如果t1是开机状态不管t1但我们假设t1停机所以t1开机就是动作则前T_min_up保持1 if t 1 Constraints [Constraints, implies(u(i,1)1, u(i,1:min(gen_params(i,6),nt))1)]; end end % 最小停机时间同理 for t 1:nt if t gen_params(i,7) - 1 nt continue; end Constraints [Constraints, implies(u(i,t-1)-u(i,t)1, u(i,t:tgen_params(i,7)-1)0)]; % 补t1的最小停机约束如果t1是停机状态同样假设t1停机其实这里t1停机的话动作不存在所以可以不用补但如果初始状态不一样的话要调 end end % 求解 ops sdpsettings(solver,cplex,verbose,1); result optimize(Constraints, cost, ops); if result.problem 0 disp([Case2 最优总成本, num2str(value(cost)), 元]); % 对比一下case1和case2的成本肯定case2高因为多了硬约束嘛 disp(Case2 机组状态矩阵); disp(value(u)); else disp(Case2 求解失败); end刚才的代码里t1的时候处理了u(i,0)不存在的问题——直接默认t1的状态就是初始动作触发的结果用implies(u(i,1)1,...)就行因为如果t1是开的那这个约束就太严了但新手入门假设初始全停是没问题的之后想改的话可以加个初始状态变量u0(i)然后把所有t1的前后差约束改成和u0(i)有关的。对比case1和case2的结果你会发现case2里的机组不会像case1那样跳来跳去比如小机组可能不会在负荷低的时候频繁启停大机组可能开机后要一直撑到满足最小时间才能停总成本也会涨个几千块这就是硬约束带来的“代价”。接下来是case3加量化功率平衡不对等下原摘要说的是“量化功率平衡”哦可能是翻译或者笔误应该是“旋转备用约束”吧电力系统里必须留一部分旋转备用容量防止负荷突然涨了或者某台机组突然跳了保证系统安全。旋转备用一般分正负不过入门练手可以只加正旋转备用所有开机机组的Pmax - 实际出力之和要大于等于某个值比如取当时负荷的10%或者加上最大单台机组的Pmax后者更安全N-1原则的旋转备用简化版。39-机组组合调度代码yalmipcplex 摘要代码是一个用来学习机组组合、熟悉yalmipcplex编程的好案例按照机组组合的约束以及难度依次设置了四个不同的case案例1为具有线性开关的简单机组组合案例2为在案例1的基础上添加了最小开关机时间约束案例3在2的基础上加入了量化功率平衡案例四则在案例3基础上实现了多天多次循环模拟 4、实现效果见下图对原摘要可能是输入错了“量化”听起来怪怪的旋转备用是UC里非常典型的“量化约束”——要留够多少量的备用所以我们就按这个来写顺便提一句原摘要可能的笔误。% --- Case3: Case2基础上加N-1简化版正旋转备用约束 --- yalmip(clear); P sdpvar(ng, nt, full); u binvar(ng, nt, full); cost 0; % 目标函数照搬 for t 1:nt for i 1:ng cost cost gen_params(i,5)*P(i,t); if t 1 cost cost gen_params(i,3)*max(0, u(i,t)-u(i,t-1)); cost cost gen_params(i,4)*max(0, u(i,t-1)-u(i,t)); end end end Constraints []; % 功率平衡和线性开关、最小开关机时间照搬case2 for t 1:nt Constraints [Constraints, sum(P(:,t)) load(t)]; for i 1:ng Constraints [Constraints, P(i,t) gen_params(i,1)*u(i,t)]; Constraints [Constraints, P(i,t) gen_params(i,2)*u(i,t)]; end end % 最小开关机时间约束和case2完全一样这里为了省篇幅就复制粘贴过来实际写代码可以写成函数四个case调用就行新手看着顺就直接放 for i 1:ng for t 1:nt if t gen_params(i,6) - 1 nt continue; end Constraints [Constraints, implies(u(i,t)-u(i,t-1)1, u(i,t:tgen_params(i,6)-1)1)]; if t 1 Constraints [Constraints, implies(u(i,1)1, u(i,1:min(gen_params(i,6),nt))1)]; end end for t 1:nt if t gen_params(i,7) - 1 nt continue; end Constraints [Constraints, implies(u(i,t-1)-u(i,t)1, u(i,t:tgen_params(i,7)-1)0)]; end end % --- 新增N-1简化版正旋转备用约束 --- % 旋转备用要求sum(u(i,t)*(Pmax(i)-P(i,t))) max(Pmax(i)) 0.05*load(t) % 这里取最大单台机组的最大出力防止这台机跳了 负荷的5%防止负荷超预期新手可以自己调比例 reserve_req (t) max(gen_params(:,2)) 0.05*load(t); for t 1:nt Constraints [Constraints, sum(u(:,t).*gen_params(:,2) - P(:,t)) reserve_req(t)]; end % 求解 ops sdpsettings(solver,cplex,verbose,1); result optimize(Constraints, cost, ops); if result.problem 0 disp([Case3 最优总成本, num2str(value(cost)), 元]); disp(Case3 机组状态矩阵); disp(value(u)); % 可以顺便输出一下每个小时的旋转备用是否够验证约束 disp(Case3 各小时实际旋转备用); disp(sum(value(u).*gen_params(:,2) - value(P), 1)); else disp(Case3 求解失败); end验证约束的时候新手一定要做比如这里输出各小时的实际备用你会发现它一定大于等于我们设的reserve_req(t)这说明约束加对了。加了备用之后总成本会再涨因为系统需要留更多的“空闲”容量可能会多开一台小机组或者让大机组的出力离Pmax远一点。最后是case4多天多次循环模拟。这个很实用比如要模拟一周的机组组合每天的负荷曲线可能是一样的典型日也可能不一样比如工作日和周末这里我们先假设每天的负荷曲线都一样就是开头给的load循环3天原参数里的case4_days。多天模拟的核心是什么是前一天的结束状态要作为后一天的初始状态不能再像前三个case那样每天初始全停了不然周一最后一台开着的大机组周二早上突然全停再重新开机成本会爆炸也不符合实际。那怎么把前一天的最后一个小时的状态u(i,nt)和后一天的第一个小时的状态u(i,nt1)联系起来还有最小开关机时间约束在两天交界的地方也要生效比如周日晚上22点开机的机组最小开机时间是4小时那周一凌晨0点、1点它也必须开着不能因为跨天了就不管了。所以这里我们的时间变量要扩展成totalnt case4days * nt负荷也要扩展成重复case4days次的load。最小开关机时间约束的循环范围也要覆盖到totalnt不用单独处理交界的地方因为循环是连续的自然就把交界的状态考虑进去了还有第一天的初始状态怎么办新手可以假设第一天之前的Tminup-1和Tmindown-1小时内所有机组都是停机状态或者更严谨一点假设第一天第一个小时之前的状态u0(i)等于第一天最后一个小时的状态u(i,total_nt)——也就是周期性边界条件这样模拟多天的循环结果会更稳定不会因为第一天的初始状态太随意而影响整体。这个周期性边界条件很有意思我们可以加上试试。% --- Case4: Case3基础上实现3天多次循环模拟 周期性边界条件 --- yalmip(clear); total_nt case4_days * nt; % 总小时数72 load_total repmat(load, 1, case4_days); % 重复3天的负荷 P sdpvar(ng, total_nt, full); u binvar(ng, total_nt, full); cost 0; % 目标函数同样覆盖所有72小时启停成本也要考虑跨天的比如第24小时和第25小时之间 for t 1:total_nt for i 1:ng cost cost gen_params(i,5)*P(i,t); if t 1 cost cost gen_params(i,3)*max(0, u(i,t)-u(i,t-1)); cost cost gen_params(i,4)*max(0, u(i,t-1)-u(i,t)); end end end Constraints []; % 功率平衡和线性开关约束覆盖total_nt for t 1:total_nt Constraints [Constraints, sum(P(:,t)) load_total(t)]; for i 1:ng Constraints [Constraints, P(i,t) gen_params(i,1)*u(i,t)]; Constraints [Constraints, P(i,t) gen_params(i,2)*u(i,t)]; end end % 最小开关机时间约束覆盖total_nt连续循环自然处理跨天 for i 1:ng for t 1:total_nt if t gen_params(i,6) - 1 total_nt continue; end % 这里不用单独处理t1了吗不对因为我们有周期性边界条件 % 周期性边界条件u(i,0) u(i,total_nt)u(i,total_nt1) u(i,1) % 所以t1的前后差是u(i,1)-u(i,total_nt)ttotal_nt的前后差是u(i,total_nt)-u(i,total_nt-1) % 那implies里的t1的u(i,t-1)就是u(i,total_nt)怎么表示 % 哦对可以先声明一个“伪”前序变量u_prev(i)然后让u_prev(i) u(i,total_nt)再把t1的约束单独拿出来用u_prev % 或者直接把implies里的u(i,t-1)在t1的时候替换成u(i,total_nt)用if-else判断 if t 1 Constraints [Constraints, implies(u(i,t)-u(i,total_nt)1, u(i,t:tgen_params(i,6)-1)1)]; else Constraints [Constraints, implies(u(i,t)-u(i,t-1)1, u(i,t:tgen_params(i,6)-1)1)]; end end % 最小停机时间同理 for t 1:total_nt if t gen_params(i,7) - 1 total_nt continue; end if t 1 Constraints [Constraints, implies(u(i,total_nt)-u(i,t)1, u(i,t:tgen_params(i,7)-1)0)]; else Constraints [Constraints, implies(u(i,t-1)-u(i,t)1, u(i,t:tgen_params(i,7)-1)0)]; end end end % --- 新增周期性边界条件可选但推荐 --- % 其实刚才的最小开关机时间约束已经用到了u(i,total_nt)作为u(i,0)但我们可以再显式加一个不用刚才的implies已经绑定了 % 不过如果想让第一天和最后一天的状态更连贯比如出力也尽量差不多不过成本函数会自然优化的显式加不加都行 % N-1简化版正旋转备用约束覆盖total_nt reserve_req_total (t) max(gen_params(:,2)) 0.05*load_total(t); for t 1:total_nt Constraints [Constraints, sum(u(:,t).*gen_params(:,2) - P(:,t)) reserve_req_total(t)]; end % 求解这里总变量数变多了72*3216个连续变量216个二进制变量不过CPLEX还是能秒算 ops sdpsettings(solver,cplex,verbose,1); result optimize(Constraints, cost, ops); if result.problem 0 disp([Case4 3天总最优成本, num2str(value(cost)), 元]); disp([Case4 平均每天最优成本, num2str(value(cost)/case4_days), 元]); % 可以画个图看看3天的负荷、各机组出力和状态不过原摘要说“实现效果见下图”这里就不放画图代码了新手可以自己用plot加subplot画 disp(Case4 最后一天的机组状态矩阵验证周期性边界条件看看和第一天是不是差不多); disp(value(u(:, (case4_days-1)*nt1 : case4_days*nt))); else disp(Case4 求解失败); end刚才的周期性边界条件真的很妙——不用单独设置第一天的初始状态直接让第一天的“前一天”等于最后一天这样求解出来的第一天和最后一天的状态会非常连贯平均每天的成本也会比每天单独初始化全停要低很多更符合实际的长期调度。总结一下这个分阶段的代码从最简单的线性开关出力约束到加最小开关机时间的硬约束再加旋转备用的安全约束最后到多天循环的实用场景每一步都只加一两个核心约束新手可以一步步敲一步步看结果的变化理解每个约束的作用同时也能熟悉YALMIP声明变量、搭约束、配求解器的流程。等下原摘要说的是“量化功率平衡”会不会我理解错了比如量化成整数兆瓦比如P(i,t)必须是整数如果是那样的话把P的sdpvar声明从full改成int就行或者int加full对sdpvar(ng, nt, int, full)。这样的话就是混合整数线性规划里的“纯整数出力”或者“整数步长出力”也算是一种“量化”。不过不管是旋转备用还是整数出力核心的分阶段思路都是一样的。