ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

港口泊位与能源系统联合优化:MILP+Gurobi实战解析

港口泊位与能源系统联合优化:MILP+Gurobi实战解析

做港口能源这行的朋友应该都有同感:港口看似是个用电大户,但真正的省钱空间往往不在设备本身,而在调度逻辑里。我做的这个项目,核心是把船舶的泊位分配和港口的综合能源系统放进同一个数学模型里联合优化——船几点到、停哪个泊位、要不要接岸电,这些决策直接决定了港口负荷曲线的形状;而光伏、储能、电网购电怎么配合,又完全取决于负荷长什么样。两边一起算,才能真正把运营成本压下来。

建模范式上我选了混合整数线性规划(MILP),求解器用Gurobi,数据处理和结果可视化用Matlab。这套组合做港口能源调度算是比较贴近工程实践的方案,也踩过不少坑。这篇文章把模型的拆解思路、数学化过程、求解器使用细节以及Matlab仿真的完整链路都整理出来,适合正在做港口减碳、综合能源调度、或者刚准备把运筹优化落地到实际项目的朋友参考。

1. 项目整体思路:为什么把泊位调度和能源系统放进同一个模型

1.1 港口负荷的“最后一公里”堵点在哪

港口电气化程度这些年提升很快,岸桥、场桥、冷藏集装箱、岸电、照明、办公负荷叠加在一起,一个中型港口的峰值动辄几十兆瓦。传统能源管理方案一般把负荷曲线当成给定输入,照着这个曲线去优化机组出力、储能充放电,最多再做个需求响应。但真正做过项目就会意识到,港口负荷根本不是“天生”的,它是调度出来的——船舶提前或延后一小时靠泊,岸电负荷就整体平移一小时;哪条船接岸电、哪条船烧柴油,直接改变电力需求的大小。

换句话说,港口负荷自带时间灵活性,而且这种灵活性的开关掌握在泊位调度手里。如果能源调度系统只被动接住一个“既定”的负荷曲线,那储能削峰填谷的空间就非常有限。反过来,泊位调度如果不看电价曲线,排出的靠泊计划很可能刚好撞上电价尖峰,买了最贵的电。这就是我在做这个项目时最先想通的一点:泊位和能量必须放在一个模型里联动,否则永远是在打补丁。

1.2 联合建模到底能带来什么收益

分开做这套问题的常规流程是:第一步排泊位,目标可能是船舶等待时间最短或者作业完成时间最早;第二步把排好的泊位计划转成岸电负荷曲线,再去做储能充放电和购电优化。逻辑上好像没问题,但两个子目标之间是断开的。

做过对比实验之后,我对联合优化的认识变得具体了很多。给一个简单的例子:某条船如果上午十点靠泊,正值光伏大发,船又能接岸电,那港口的自发电就被就地消纳了,电网购电量下降;但如果只追求泊位利用率最高,把这条船排到下午两点靠泊,正好赶上电价峰段,储能放完电还不够,只能高价买电。两种方案单独看都“合理”,合在一起看成本差一大截。联合建模通过共享决策变量,迫使模型在所有物理约束里找一个整体最优解,泊位时机、充放策略、购电曲线全部联动起来。这才是“协同”两个字的本意。

1.3 技术选型:为什么是MILP加Gurobi,Matlab又扮演什么角色

港口泊位分配天然包含离散决策,比如某条船在某个时刻是否占用某个泊位、是否接入岸电、储能是否处于充电状态,这些都得用二进制变量表达。而能量流、功率、SOC都是连续量,可以用线性方程描述。两个部分拼在一起,问题的标准形态就是混合整数线性规划,也就是MILP。

选Gurobi首先是看中它的求解性能和稳定性。做这类工程优化最烦的就是模型稍微大一点就求解半天,或者数值稍有问题就疯狂报错。Gurobi在MIP求解器里属于第一梯队,自带的不可行性分析IIS功能调试约束冲突也特别好用。Matlab在这套流程里主要承担数据前处理和结果可视化。我习惯用Matlab做仿真算例的搭建、画泊位甘特图、画功率堆叠图,比Python的matplotlib省事不少;核心优化求解则交给Gurobi,两边通过数据文件交换信息。这样分工最大的好处是:哪一侧出问题都好排查。

2. 泊位分配子模型:把“船靠哪”变成数学语言

2.1 决策变量与参数怎么设计

泊位分配子模型的第一步是时间离散化。我把一个调度周期分成T个时槽,最常见的做法是1小时一个点,需要更精细控制时可以细化到15分钟,但代价是变量数量成倍增长。设参与调度的船舶集合为V,泊位集合为B,定义几个关键决策变量:

  • z[v][b][t]:二进制变量,1表示船v在t时槽占用泊位b;
  • a[v]:船v的实际靠泊时刻;
  • d[v]:船v的实际离泊时刻;
  • y[v]:二进制变量,表示船v是否在靠泊期间接入岸电。

输入参数方面,每条船带一个预计到港时间窗[E[v], L[v]],还有一个作业时长p[v]。这里要注意,作业时长不只是装卸货的时间,还包括抛缆、接岸电、安全检查这些辅助流程。泊位这边还要区分类型,有的泊位带高压岸电、有的带低压岸电、有的干脆没有岸电设施,这些都要作为参数拍进模型。

2.2 核心约束的构造逻辑

扎约束之前先说一个我踩过的坑:如果用z[v][b][t]来表达“一条船占用某个泊位且连续作业”,最直观的写法是sum_t z[v][b][t] == p[v],但这样只能保证总时长为p[v],并不能保证时段连续。船在泊位上干两小时歇一小时再干两小时,这种物理上不可能的方案,在这种写法下反而是可行解。

我最终的解决办法是引入靠泊时刻和离泊时刻,把连续占用的判断交给人这个模型去推导。具体做法是:如果a[v] <= t <= d[v],那么z[v][b][t]必须等于某个泊位b的占用标识。变量之间的绑定用大M约束实现,比如z[v][b][t] >= 1 - M * |t - a[v]|这种写法。虽然严谨的等价变换需要几个约束组合,但核心思想就是强制让占用区间和靠离泊时刻对齐。

除此之外,最基础的约束还有三条:每条船只能占用一个泊位;任何一个时刻同一个泊位最多只能有一条船;每条船的靠泊时刻必须在到港时间窗内。第三条约束在实际工程中往往要加松弛变量,因为船舶到港受天气、航道、前一港拖班影响很大,硬性约束会让模型经常无解。我在项目里给时间窗加了一个可惩罚的松弛变量,也就是允许晚靠,但晚靠多少会在目标函数里扣分,这样模型既有可行性,又能体现优先级。

2.3 岸电负荷如何“长”在泊位上

岸电负荷是整个港口综合能源系统里最核心的耦合量。它由泊位分配直接决定,同时又作为能源子模型的输入参与电力平衡。我在模型里这样写:

P_shore[t] = sum_v sum_b z[v][b][t] * y[v] * P_need[v]

其中P_need[v]是船v在接岸电状态下的电力需求。这个公式看起来很直白,但要注意一个工程上的隐形细节:船接了岸电之后,并不是靠泊全过程都在满功率用电,它有靠港初期的启动峰值、作业中的稳定负荷、离港前的逐步退出。我在模型里用的是平均功率近似,如果要做高精度仿真,可以把每个船的用电曲线按时间分段,作为已知序列和一串二进制绑定,问题规模会大幅上升,但思路完全一致。

还有一个细节:接岸电是船方的选择,不是强制要求。实际运营中,接岸电往往伴随碳排放配额收益、船舶服务费优惠,这些可以量化为目标函数里的线性项。把这些收益写进优化目标之后,模型才会主动“引导”船接岸电,而不是因为它只是约束里的一个可选项就被丢弃。

3. 多能协同能源系统子模型:光储充怎么配合

3.1 设备数学模型:光伏、储能、电网三个模块

能源系统这边我最关心的设备是光伏、储能和电网接口,它们构成港口电力母线的三个主要电源。

光伏模型最简单直接:给定预测出力序列P_pv[t],如果允许弃光,就加一个弃光变量P_cut[t],实际出力等于预测值减弃光量。不要小看弃光变量,很多场景里光伏大发遇到负荷低谷,储能容量又不足以吸收全部电量,这时候允许弃光能让模型有解,也反映真实运行中的限电操作。

储能模型的重点是SOC递推方程。我用的经典写法是:

E[t] = E[t-1] + (P_ch[t] * eta_ch - P_dis[t] / eta_dis) * dt

其中P_ch和P_dis分别是充放电功率,eta_ch、eta_dis是充放电效率。这里有一个新手容易犯错的地方:直接在等式里让储能同时充电又放电,数学上模型是允许的,因为充放抵消表面上不亏能量,但实际中同时充放会白白产生损耗,对运行成本非常不友好。我加了互斥约束,用两个二进制变量分别控制充放状态,让它们不能同时为1。另外SOC的初始值和末端值也要给约束,比如要求最后一个时段的SOC不得低于初始值的50%,否则模型会为了贪小便宜把储能彻底放空,仿真结果毫无可持续性。

电网购电模块相对简单,功率上下限对应变压器的容量约束即可。如果厂区里还接入了电制冷、制热设备或者氢能设备,建模方式类似,本质都是在能量平衡方程里增加一项负荷或一个转换单元。

3.2 能量平衡与耦合约束是两个模型的黏合剂

把泊位子模型和能源子模型缝合到一起的,就是一条覆盖所有时段的功率平衡约束:

P_grid[t] + P_pv_actual[t] + P_dis[t] = P_base[t] + P_shore[t] + P_ch[t]

等号左边是电网购电、光伏实际出力和储能放电,右边是港口固定负荷、岸电负荷和储能充电功率。只要这条约束成立,整条母线上的功率就守恒了。P_shore[t]正是从泊位子模型里导出的岸电功率,它是这两个子模型发生交互的唯一桥梁,也是整个联合优化模型最关键的一条等式。

这条约束写起来只有一行,但实际调试时最花时间。因为任何一个子模型里变量的单位不一致、时间索引偏移或者功率正负方向约定不一,最后都会反映到这条平衡约束上。比如储能放电我定义成正数,但有人习惯定义成负数,一旦混用,平衡约束就乱套。建议在建模初期就统一所有变量的正方向约定,我通常把所有输出功率定义为正方向。

3.3 目标函数的设计取舍:成本还是碳排

目标函数我选了最小化综合运行成本,组成分三块:电网购电费用、船舶等待时间惩罚、储能充放电损耗近似惩罚。购电费用按分时电价累加,分时电价越高,模型越趋向于错峰用电,这是整个调度的经济指挥棒。船舶等待时间惩罚用加权系数表示优先级,重要客户船或冷冻货船权重高,晚靠的代价就大,模型会在泊位分配时优先照顾它。储能充放电损耗惩罚是一个很小的线性项,目的是避免无意义的充放循环。

为什么不直接最小化碳排放?我做过对比实验,目前多数港口运营方考核的硬指标还是钱,碳排放更多是以约束或权重形式出现。如果项目要求优化碳排,备选方案是把目标函数改成最小化碳排放,同时把购电成本写成约束上限,由决策者给定一个可接受的成本预算。两种写法各有优缺点:加权重简单直观但对系数敏感,分层约束物理意义更清晰但可能牺牲一部分经济性,实际项目里需要跟业主确认到底哪个指标更硬。

4. Gurobi求解落地:从数学模型到能跑的代码

4.1 用gurobipy搭建求解骨架

Gurobi的原生Python接口gurobipy在建模上比Matlab的YALMIP更直接,而且调试报错信息非常清楚。我通常把整个求解代码分成数据读入、模型构建、求解、结果导出四个模块,下面是模型构建部分的核心骨架:

import gurobipy as gp from gurobipy import GRB m = gp.Model("port_energy_coopt") # 泊位分配变量 z = m.addVars(V, B, T, vtype=GRB.BINARY, name="z") y = m.addVars(V, vtype=GRB.BINARY, name="y") a = m.addVars(V, vtype=GRB.INTEGER, lb=0, ub=Horizon, name="arrive") d = m.addVars(V, vtype=GRB.INTEGER, lb=0, ub=Horizon, name="depart") # 储能与电网变量 e_soc = m.addVars(T + 1, lb=0, ub=E_max, name="soc") p_ch = m.addVars(T, lb=0, ub=PB_max, name="charge") p_dis = m.addVars(T, lb=0, ub=PB_max, name="discharge") u_ch = m.addVars(T, vtype=GRB.BINARY, name="u_ch") u_dis = m.addVars(T, vtype=GRB.BINARY, name="u_dis") p_grid = m.addVars(T, lb=0, ub=P_grid_max, name="pgrid") p_cut = m.addVars(T, lb=0, ub=P_pv_max, name="curtail") p_shore = m.addVars(T, lb=0, name="pshore")

约束的写法这里不全部展开,就说一个我常用的技巧。当模型里既有连续变量又有大量二进制变量时,条件约束的表达式非常繁琐。比如“如果船v不在泊位b上,那么该泊位对应的岸电功率为0”这条约束,我惯用的做法是把所有逻辑判断都转化成大M不等式。大M的值能取多小就取多小,这直接决定求解器数值稳定性。比如时间量纲是小时,M取24就够了,千万别无脑给1E6,否则既拖慢求解又容易报numerical issues。

4.2 求解参数整定:让模型在可接受时间内给出好解

Gurobi就算法本身非常强,但参数不会调一样可能几个小时跑不出结果。我一套比较稳的配置是这样:

m.Params.MIPGap = 0.005 m.Params.TimeLimit = 900 m.Params.Threads = 4

MIPGap是求解器停止的精度门槛,0.005意味着最优解和当前解的相对差距小于千分之五就收手。工程场景下这个精度完全够用,非要追求0.0001只会白白消耗算力。TimeLimit设900秒,也就是15分钟,到点就把当前最优可行解吐出来。实际接入调度系统时,这个时间必须卡得住,否则没法做滚动优化。

数值稳定性是我在这个项目里花时间最多的地方。第一次用kW、kWh、元/kWh混合建模,目标函数值很大,Gurobi的容差默认值相对就变小了,判断整点、比较上下界的逻辑都容易出问题。后来把所有物理量统一成MW和MWh,目标函数数值落在一个温和的区间,MIP的收敛速度明显改善。这是做这类项目的一条铁律:单位必须全局统一,最好在代码入口处一次性换算。

4.3 结果导出与Matlab联动

求解完成之后拿到变量值,我习惯用m.getAttr('X', z)把二进制变量的矩阵取出来转成numpy数组,再通过scipy.io.savemat存成Matlab能读的.mat文件。这样Matlab侧直接load就能拿到泊位占用矩阵和储能SOC序列,作图特别顺手。

如果团队更习惯全Matlab流程,用YALMIP工具箱包一层也可以,核心调用只要几行:

ops = sdpsettings('solver', 'gurobi', 'verbose', 2); optimize(Constraints, Objective, ops);

但这里有个工程判断:YALMIP作为建模层很方便,可一旦模型规模变大,中间层引入的变量和约束膨胀也比较明显。我自己的经验模型规模大时还是直接用gurobipy更可控,Matlab专心做前处理和可视化,分工清晰,调试效率最高。

5. Matlab仿真与算例分析:数字到底说明了什么

5.1 仿真场景与数据准备

为了验证模型,我搭了一个典型的中型港口算例。仿真周期24小时,时间步长1小时,一共3个泊位、5条船。船的参数如下表,这里面的数据是示意用的非实测值,但形状和量级来自真实调度报表的压缩。

船编号到港窗口作业时长泊位类型需求岸电功率
V10-4点4小时高压岸电800 kW
V28-14点2小时低压岸电400 kW
V310-20点6小时不需要0
V412-22点3小时低压岸电500 kW
V520-24点5小时高压岸电600 kW

分时电价采用典型的峰谷结构:峰段(10-15点、18-21点)1.2元/kWh,平段0.7元/kWh,谷段(0-6点)0.35元/kWh。光伏功率曲线用一条晴好天气的典型日曲线,峰值出现在正午,装机容量300 kW。储能规格是600 kWh容量、300 kW功率,充放电效率取0.95。

数据准备阶段,我强烈建议在Matlab里做一张数据字典表,写清楚每个序列的时间索引是绝对时刻还是相对时槽,单位是kW还是MW。这个项目里我在这上面吃过亏——有一次Matlab画出的负荷曲线和Python侧算出来的完全对不上,排查了一个多小时才发现是一个数据文件里混用了两种时间基准。

5.2 泊位甘特图与能量曲线怎么读才有价值

仿真结果出来之后,有两张图是我每轮必看的:泊位甘特图和源-荷功率堆叠图。

泊位甘特图的横轴是24小时,纵轴是三个泊位,每个色块代表一条船在某个泊位的占用区间。这张图能直接看出模型调度的“心思”。我在算例里观察到的一个很有意思的现象:V5这条船的窗口很晚,正常情况下会排到夜间靠泊,但联合优化后的结果里,它被安排到凌晨时段靠泊,刚好赶在低谷电价时段接岸电。表面上看,凌晨靠泊增加了船舶等待成本,但对比夜间谷电和傍晚峰电的价差,综合成本反而更低。这种“牺牲一点泊位效率、换来整体费用下降”的决策,正是协同优化最直观的价值。

功率堆叠图则是把P_pv_actual、P_dis、P_grid画成堆叠柱状,把P_base、P_shore、P_ch画成负荷堆叠。两相对照,能很清楚地看到储能什么时候充电、什么时候放电,光伏出力去了哪里。我重点看SOC曲线的走势,好的模型SOC曲线应该呈现低电价充电、高电价放电的锯齿形,且放电尖峰正好顶在岸电负荷最高的时候。如果SOC曲线和电价走势没有同步关系,说明模型约束哪里出了问题,白白浪费了储能的灵活性。

5.3 联合优化对比分步优化:省下来的钱从哪来

我把算例拆成两种方案做了一组对照实验。方案A是传统分步流程,先最小化船舶等待时间得到泊位计划,再把这个计划固定下来做能量调度。方案B是本文的联合优化。结果挺直观:方案B的综合运行成本比方案A下降了约9.6%。

成本的下降主要来自三个方面,我拆开算过。第一,泊位分配避开了高电价时段,V5和V1两条大功率岸电船错峰靠泊,光这一项就省下了可观的购电费用。第二,储能充电时机对齐了光伏出力,午间光伏大发时储能充满,傍晚高峰时放电,替代了高价电网电。第三,岸电接入决策本身参与了优化——有的船在峰电时段干脆不接岸电,优先保证更紧要的作业用电,这个策略在分步流程里根本没法实现。

这个9.6%是在一个小型算例上跑出来的,如果放到实际港口,船舶数量更多、负荷波动更大,错峰与协同的空间通常只会更大。当然,联合优化模型规模大,求解时间也长,实际部署时要结合调度时间窗口做取舍。

5.4 模型扩展方向:从电力母线到多能联供

做完这套基础模型,扩展空间其实很大。港口不只是用电,还有冷库制冷、办公采暖、甚至制氢加氢的用氢需求。把电力母线的模型扩展到电-热-冷联供,只需要在能量平衡方程里增加热母线、冷母线的对应等式,新增电锅炉、吸收式制冷机、储热罐这类设备的线性化模型,然后让它们共享泊位分配的决策信息。我实际试过在模型里增加一个电制冷机模块,增加的求解时间不大,但能量梯级利用带来的成本下降很可观。

另一个有价值的扩展方向是滚动时域调度。把24小时拆成若干个窗口,每4到6小时重算一次,保留最后几个时段的决策作为实际执行指令。这种方式特别适合应对船舶晚点、天气突变这些实时扰动,是模型从离线算例走向现场调度的必经一步。

6. 常见问题与排查技巧实录

6.1 模型报无解:先查IIS,别靠猜

做这套模型,最先遇到的基本都是无解。没有经验的人看到infeasible状态就慌了,开始随机删约束。我建议第一步先跑m.computeIIS(),再用m.write("model.ilp")把冲突集合写出来,查看被标记的约束是哪些。

我在算例里踩过最典型的一个坑是时间窗和作业时长互相矛盾:某条船的预计到港窗口结束时刻减去作业时长,已经不够完成作业了,比如窗口是20到24点,作业时长6小时,这船根本排不进去。这种问题如果不是靠IIS一眼定位,纯靠人工审查约束根本查不过来。另一个常见问题是大M值的符号方向写反,导致某条约束永远不可能满足。遇到无解,先查IIS永远比乱试约束高效。

6.2 求解时间爆炸:三个手段依次上

模型跑不动是常态,想到几个缓解手段。第一步调MIPGap,从0.001放到0.02甚至0.05,很多工程场景百分之五以内的次优解完全够用。第二步加对称破除约束,当多个泊位能力完全相同时,模型会产生大量等价解空间,搜索效率极低。比如强制船舶按编号顺序从编号较小的泊位开始分配,能显著压缩分支数。

如果这两个手段还不够,就上滚动时域。把24小时拆成4个6小时窗口,每个窗口之间的重叠时段设置1小时,把上一轮末尾的SOC、泊位占用状态作为下一轮的初始条件,这样每个子问题规模变小,求解时间从小时级降到分钟级。代价是失去了全局最优,但工程调度里实时性往往比一点次优损失更重要。

6.3 与Matlab联调的几个常见坑

Matlab和Gurobi联调时最容易翻车的是单位不一致。比如Matlab读入的光伏预测数据是kW,Python侧模型单位是MW,平衡方程里所有的量表面上都能对上,但结果完全错乱。我在数据文件头部写清单位标注,在代码入口统一除以1000换算成MW。

另一个高发问题是时间索引偏移。Python数组从0开始,Matlab从1开始,两者交换数据时在循环里少写一个偏移,结果整条曲线错位一小时。更隐蔽的是出现在跨时段的SOC递推约束里,错位一个时槽,SOC曲线就会产生诡异的突变。我的习惯是为所有时间序列配一个显式的绝对时间轴变量,存储和读取都按绝对分钟数对齐,程序里不再依赖数组下标做时间对应。

还有一个容易被忽略的点是随机数种子。如果给储能初始SOC、船舶到港时间这些参数加了随机扰动,而没有固定随机种子,每次跑出来的结果都不一样,做对照实验时会很头疼。在Matlab用rng(2024)固定随机种子,Python侧用np.random.seed(42),所有实验都具备可复现性。

6.4 一个容易被忽略的建模细节:SOC与功率的取整

储能SOC递推方程里,E[t]和E[t-1]是连续变量,在求解器看来,任意小的充放电变化都会带来SOC的小幅调整。如果求解器数值容差没控制好,会出现SOC曲线在约束边界上反复跳动的情况,这在Matlab图里看起来就是一条锯齿状的线。

这个问题的根源是模型里变量数量过多、数值梯度太密。我的解决办法是给SOC加一个很小的正则项,逼迫模型不要做无谓的SOC抖动,同时对所有SOC变量设置合理的容差范围。这样做之后结果曲线平滑很多,和实际情况也符合得更好。这个细节属于那种不会导致模型错误、但会让仿真结果看起来不专业的问题,值得记录一下。

这套项目走下来,我最深的体会是:港口能源优化的难点真的不在算法,而在把所有物理逻辑准确翻译成数学约束的过程。泊位分配和能量调度看起来是两个专业方向,但放到同一个优化框架里,反而能产生单模块优化永远给不出的协同价值。建议想上手的朋友不要一开始就追求模型完美,先把1小时步长、3个泊位、5条船的小算例跑通,把每条约束背后的物理含义都理清楚,再逐步加碳排放约束、多能联供、滚动时域这些进阶功能。模型规模和求解时间永远是一对矛盾,但先小后大、先通后优这条路基本不会走偏。

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进