
前段时间我在做园区多能互补调度项目的收尾时被一个现象折腾得够呛用固定效率、固定碳价算出来的储能充放电策略拿到实际现场跑了一周综合成本比模型预测值高出将近8%。不是模型框架出了问题而是两个被我当成常数的东西在真实系统里根本就不是常数——一个是储能损耗一个是碳价对成本的边际影响。今天把这套改进模型完整拆开讲一遍。这个模型做的事情很简单在传统多源协同调度基础上把储能损耗从单一效率换成按功率分段的损耗曲线把碳交易成本从固定单价换成阶梯碳价又把需求侧响应作为可调资源一起塞进优化。难点不在各自单独建模而在于分段损耗、阶梯碳价、需求响应报价这三类非线性约束放进同一个混合整数规划里会互相拖累求解速度也更容易踩数值坑。如果你正在用Python写储能优化、综合能源调度、碳成本相关模型或者只是想知道怎么把这类分段线性问题干净利落地写进求解器这篇笔记应该能帮到你。1. 为什么固定碳价和恒定储能效率在工程现场站不住脚1.1 固定碳价与实际账单里的边际跳变很多论文里的碳成本处理都是总排放量乘一个固定单价画出来是一条直线。这个写法本身没有错但放到真正的碳市场结算规则里它漏掉了一个关键特征碳价往往是分阶梯的配额内是一个价格超配额以后超过越多单价越贵。我给园区做模型的时候一开始也用了固定碳价结果财务那边拿着账单来问我为什么这个月排放只多了3%碳成本却多了17%因为那个月刚好跨过了配额台阶多出来的部分被更高档的单价计费出现了特别明显的边际跳变。固定碳价模型无论如何都解释不了这个现象。所以要准确估计碳成本就得把碳价做成阶梯函数。阶梯碳价不是简单加个分段线性项就完事它本质上带有整数选择变量因为系统会“选择”自己落在哪个碳价区间。这也是它比固定碳价建模麻烦的核心原因。1.2 恒定效率与实际储能损耗的偏差储能模型里的效率常数问题类似。锂电池的充放电损耗和功率不是线性关系低功率、低SOC状态下内阻变化、热损耗占比都不同做得细一点还有功率转换装置在部分负载下效率明显下降。用单一效率常数等于把所有充放电过程都当成额定工况。在我测试的算例里同一组历史运行数据用恒定0.9效率算出来的日充电量需求和用分段损耗曲线算出来的结果相差6.5%左右。别小看这个数调度是全天24小时滚动做的6.5%的偏差经过储能SOC累积经常在夜间时段触发不了原本应该发生的放电动作最终导致负荷高峰时还要高价购电。分段损耗的思路是把充放电功率拆成几个功率区间每个区间用不同的近似效率或者用一条折线函数去逼近真实损耗曲线。这种处理能抓住“大功率充电损耗高”这个行业常识又保留了线性模型的数学性质代价只是多加点辅助变量。2. 模型主框架多源协同里每一个成本项怎么进目标函数2.1 决策变量与约束体系这个模型的时间尺度是24小时时间步长1小时。对象包括风电、光伏、燃气轮机、电储能、电网购电以及需求侧可削减负荷。目标函数是系统总成本最小包含外购电成本、燃气机组燃料成本、储能运维成本、需求响应补偿成本、碳交易成本。决策变量分成四组第一组是各电源出力第二组是储能充放电功率和SOC第三组是需求响应削减量第四组是为了处理分段损耗和阶梯碳价引入的辅助0-1变量。千万别小看最后一组变量阶梯碳价和分段损耗的整数变量加在一起会让整个模型从线性规划变成混合整数规划求解难度完全不在一个量级。很多初学者喜欢把碳成本直接丢进目标函数然后加一个功率平衡方程就完事。但真正的难点在约束之间的耦合关系功率平衡要把储能充放电、需求响应削减、电网购电都包含进去SOC约束要跟分段损耗联动碳排放约束又依赖燃气机组出力和电网购电量的乘积。这些约束是相互咬合的每一步都会影响储能要不要充电、充多少、什么时候放。2.2 分段损耗的三段式线性化处理我实现里用的是三段式损耗充电功率0到100 kW、100到200 kW、200到250 kW三段对应存储效率0.92、0.90、0.87放电对应输出效率0.93、0.91、0.88。这些数字来自电池厂商提供的一部分曲线我在实际算例里又做了小幅修正。直接把效率乘在功率变量上没问题问题出在分段选择。功率落在哪一段本身是一个离散决策因为目标函数希望压低损耗会故意把功率卡在区间边缘。这时候不能用普通连续变量硬来表达否则求解器会认为它可以在两个分段之间“插值”产生一种并不存在的中间效率优化结果就会失真。处理办法有两种一种是SOS2权重变量一种是二进制变量加Big-M。Gurobi里我倾向于SOS2因为约束数量少、数值稳定性好如果用的是CBC这类开源求解器SOS2支持较弱就得退回到二进制大M法。实际工程中先把损耗曲线拟合成3到5段就够用了没必要追求几十个段点后面我会单独说这个求解时间问题。2.3 阶梯碳价引入整数变量的建模技巧阶梯碳价我采用的参数是免费配额5000 kg/天超排量0到2000 kg按0.1元/kg2000到5000 kg按0.15元/kg5000 kg以上按0.2元/kg。如果配额没用完模型里不会产生卖配额收益这个处理更贴近多数用户侧项目的真实结算方式。阶梯碳价建模的核心是把超排量拆成三段增量每一段量乘对应单价。逻辑上依赖0-1变量来保证区间使用顺序只有在第一段填满之后第二段才可能有值。这个顺序约束经常被漏掉漏掉的结果就是求解器会把超排量全部塞到单价最低的第一段去算出来的碳成本低得离谱。我在实际代码里用的是一种更稳的写法先拆出免费配额用量和富余量再把超排量拆成三个区间变量每个区间变量有上下界并用二进制变量控制激活状态。由于碳单价递增目标函数会自动优先填满低价区间所以填满约束不是必须的但为了避免在某些极端约束组合下出现反直觉解我还是保留了顺序激活条件。3. Python实现从数据预处理到求解器配置的完整链路3.1 工具链与环境准备我习惯的环境是Python 3.9加Gurobi 10搭配pandas做数据清洗、matplotlib做结果可视化。Gurobi支持SOS2这对分段损耗建模省了很多事。如果你没有Gurobi学术许可开源求解器HiGHS或者CBC也能跑但MIP求解性能差距在规模变大后会非常明显。代码文件结构我会拆成main.py、data_loader.py、model_builder.py、post_process.py模型构建和结果可视化严格分开。这里有个血泪教训不要把所有东西堆在一个脚本里。这个模型本身不算复杂但当你同时调试分段损耗、阶梯碳价、需求响应三个特征时如果全部混在一个文件里排查问题会极其痛苦。很多朋友在配置Python环境时容易被IDE的解释器识别问题卡住其实那只影响编辑器的语法提示和运行按钮模型本身只认你命令行里那个Python环境。我建议直接用conda建一个独立环境把gurobipy、pandas、matplotlib、openpyxl装好然后在PyCharm里把Project Interpreter指过去能省掉一半的环境报错。3.2 核心代码构建MIP模型先看模型骨架我用Gurobi Python API写决策变量和基础约束。import gurobipy as gp from gurobipy import GRB T 24 dt 1.0 E_bat 1000.0 # 储能容量 kWh SOC0, SOC_min, SOC_max 0.2, 0.1, 0.9 m gp.Model(mssp) # 电源与购电 P_buy m.addVars(T, lb0, nameP_buy) P_gt m.addVars(T, lb0, ub300, nameP_gt) # 需求响应削减量 P_dr m.addVars(T, lb0, nameP_dr) # 储能 P_ch m.addVars(T, lb0, ub250, nameP_ch) P_dis m.addVars(T, lb0, ub250, nameP_dis) u_ch m.addVars(T, vtypeGRB.BINARY, nameu_ch) u_dis m.addVars(T, vtypeGRB.BINARY, nameu_dis) SOC m.addVars(T 1, lbSOC_min, ubSOC_max, nameSOC) m.addConstr(SOC[0] SOC0) # 充放电互斥 m.addConstrs(P_ch[t] 250 * u_ch[t] for t in range(T)) m.addConstrs(P_dis[t] 250 * u_dis[t] for t in range(T)) m.addConstrs(u_ch[t] u_dis[t] 1 for t in range(T)) # 功率平衡 m.addConstrs( P_w[t] P_pv[t] P_gt[t] P_dis[t] P_buy[t] P_load[t] - P_dr[t] P_ch[t] for t in range(T) )P_w和P_pv在demo里直接作为已知参数传入。这里没有卖电方向所以功率平衡右边把充电当成负荷把放电当成电源。接下来是分段损耗的SOS2实现。以充电侧为例折线点横坐标是充电功率纵坐标是实际存入电池的有效功率。charge_x [0, 100, 200, 250] charge_y [0, 92, 180, 217.5] # 每个区段效率 0.92 / 0.90 / 0.87 P_store m.addVars(T, lb0, nameP_store) w_ch m.addVars(T, len(charge_x), lb0, namew_ch) for t in range(T): m.addConstr(gp.quicksum(w_ch[t, i] for i in range(len(charge_x))) 1) m.addSOS(GRB.SOS_TYPE2, [w_ch[t, i] for i in range(len(charge_x))]) m.addConstr(P_ch[t] gp.quicksum(charge_x[i] * w_ch[t, i] for i in range(len(charge_x)))) m.addConstr(P_store[t] gp.quicksum(charge_y[i] * w_ch[t, i] for i in range(len(charge_x))))放电侧同理可以用另一组SOS2变量把P_dis映射成SOC的释放量。SOC更新就可以写成m.addConstrs( SOC[t 1] SOC[t] (P_store[t] - P_release[t]) * dt / E_bat for t in range(T) )注意单位统一功率kW时间h容量kWh。如果直接把250 kW和1000 kWh放进约束变量尺度差5倍以上Gurobi的数值容差就可能出问题。后面踩坑部分我会专门展开。阶梯碳价的增量线性化代码是这样# 碳排放总量 E_gt gp.quicksum(P_gt[t] * 0.2 for t in range(T)) E_buy gp.quicksum(P_buy[t] * 0.581 for t in range(T)) E_total E_gt E_buy E_free 5000.0 E_free_used m.addVar(lb0, nameE_free_used) E_free_surplus m.addVar(lb0, nameE_free_surplus) E_excess m.addVar(lb0, nameE_excess) m.addConstr(E_total E_free_used E_excess) m.addConstr(E_free_used E_free_surplus E_free) # 超排量分三段 carbon_len [2000, 3000, 100000] carbon_price [0.10, 0.15, 0.20] z_carbon m.addVars(3, lb0, namez_carbon) y_carbon m.addVars(3, vtypeGRB.BINARY, namey_carbon) m.addConstr(E_excess z_carbon[0] z_carbon[1] z_carbon[2]) for i in range(3): m.addConstr(z_carbon[i] carbon_len[i] * y_carbon[i]) m.addConstr(y_carbon[0] y_carbon[1]) m.addConstr(y_carbon[1] y_carbon[2]) carbon_cost gp.quicksum(carbon_price[i] * z_carbon[i] for i in range(3))这段代码的精髓在于用E_free_used和E_free_surplus把“配额内”和“超配额”分开避免直接对E_total减E_free取max这种非线性操作。E_total小于免费配额时E_free_used会等于实际排放量E_excess被压到0目标函数不会为此多付一分钱。需求响应补偿可以按线性价格直接写如果补偿报价是二次函数就需要分段线性化。我建议第一版先用线性补偿把模型跑通后再考虑二次项否则排查问题时分不清是碳价阶梯的问题还是DR二次项的问题。3.3 数据输入与结果可视化数据输入我习惯用pandas读Excel结构尽量是长表每个时刻一行包含负荷、风电、光伏、电价、负荷可削减上限等。别在Excel里做复杂公式那些公式在pandas读进来之后全部失效只保留原始数值是最好的。结果可视化部分我会画三张图第一张是功率平衡堆叠图能看到每个时刻风电、光伏、燃气、购电、储能充放电的构成第二张是SOC曲线重点观察有没有频繁的充放电切换第三张是碳成本和碳排放柱状图叠加碳价台阶线。这三张图能快速定位问题如果SOC曲线一天切换超过七八次说明目标函数里少了储能寿命相关惩罚项或者分段损耗曲线没起到诱导作用。4. 算例拆解需求侧响应和阶梯碳价是怎么互相放大的4.1 四个对比场景