简介面向水风光互补调度与投资评估的Matlab代码资源适用于能源、电气工程、计算机及数学等专业学生的课程设计、期末大作业或本科毕业设计。资源以不同网格装机容量和出力系数为关键变量完整实现水风光互补系统的优化调度建模与净现值计算内容涵盖风力、光伏、水电出力估算、调度决策及经济性对比可帮助读者掌握新能源并网环境下投资效益分析的基本流程。压缩包共18个文件、约475KB其中12个M脚本承担核心算法与主流程3个Excel表格提供可直接运行的案例数据另有Python脚本、Word说明和Markdown文档辅助理解代码注释详细、参数化程度高便于修改关键参数后重新求解并观察结果变化。目前已有34人学习适合初步接触新能源调度建模、需要快速复现算例并进一步扩展实验的学生。1. 水风光互补调度里的净现值与网格分析经济账和空间尺度必须一起算拿到水风光互补调度、净现值计算、不同网格装机容量和出力系数这个.zip包时多数人的第一反应是先把调度模型跑通再把NPV算出来。但真正做过的人都知道这两件事不能分开做调度结果决定年发电量年发电量是净现值现金流的唯一收入来源而装机容量和出力系数又直接决定调度模型里各电站的出力上限和实际发电水平。网格选得粗出力系数的高估能在NPV上放大几个百分点直接把一个勉强能投的项目算成优质标的。下面我从调度模型、净现值落地、网格敏感度分析三个层面拆开讲最后列出五个我实际踩过的坑。适合正在做流域水风光一体化规划、可研阶段经济评价的工程师和数据研究人员。2. 水风光互补调度模型先把物理约束写成可求解的优化问题2.1 互补调度的核心逻辑径流、风速、辐照三条时序怎么耦合水风光互补调度本质上是一个多能源联合优化问题。水电的可调度性强水库能蓄能放风光出力受天气支配波动性大。互补的价值在于水电在风光不足时补位在风光大发时蓄水把径流资源的“时间平移”能力和风光资源的“空间散布”能力组合起来。模型里最基本的三个组件是水电的入库径流、库容水位、发电流量限制、水头-出力关系风电的切入/切出风速、额定功率、功率曲线光伏的组件容量、综合效率、辐照度到电功率的折算。调度的目标函数可以是最小化弃电率、最大化系统发电量也可以是在给定分时电价曲线下最大化收益。如果后续要做净现值我建议目标函数直接设为“最大化运营期内逐年收益的折现值”这样调度结果和财务模型做到同一个目标口径避免两套模型打架。约束条件里最容易出错的是水量平衡方程水库调度常用的一阶近似是V(t1) V(t) (Q_in(t) - Q_out(t) - Q_spill(t)) × dt其中V是库容Q_in是入库流量Q_out是发电流量Q_spill是弃水流量。弃水不是模型里的自由变量当库容达到上限时溢出的水才算弃水这需要用一个松弛变量表达新手最容易在这里翻车第5章会详细展开。风光功率的时序约束看起来简单数据质量问题往往更大。风速-功率折算建议用实际机组的功率曲线做分段线性插值而不是用一个三次方公式通吃。光伏则要把辐照度乘以组件效率、温升系数、灰尘遮挡系数、逆变器效率四个系数叠下来实际出力大约是理论值的75%到85%。这个“发电综合效率系数”的取值会直接传导到NPV的现金流里不能拍脑袋。2.2 用Python搭建最小调度模型的完整代码与参数说明下面给出一个最小可运行的日尺度互补调度求解代码用pandas做时序处理scipy.optimize.linprog做线性规划求解。为控制篇幅水电用简化模型忽略水头随库容变化风光用固定折算系数但变量排列、约束写法、单位换算是完整的可以直接替换成实际数据。import numpy as np import pandas as pd from scipy.optimize import linprog # 24小时步长水风光三站联合调度决策变量只针对水电 T 24 # 流量单位换算1 m3/s 持续1小时 0.0036 百万m3 dt 3600 / 1e6 # 入库流量(m3/s)模拟一个日调节水库的来水过程 flow_in np.array([50,48,45,42,40,38,36,35,33,32,30,29, 28,27,25,24,23,22,20,18,16,14,12,10]) # 水电参数 V_0 80.0 # 初始库容百万m3 V_min, V_max 20.0, 120.0 # 库容上下限 Q_min, Q_max 10.0, 60.0 # 发电流量上下限m3/s head 50.0 # 平均水头m eta_h 0.85 # 水轮机效率 # 出力系数rho*g*head*eta/1e6单位 MW/(m3/s) coef_h 1000 * 9.81 * head * eta_h / 1e6 # 决策变量排列[Q_out_0..Q_out_23, V_0..V_23]共48个 nvars T * 2 c np.zeros(nvars) c[0:T] -coef_h # 最大化水电出力转成最小化负值 # 不等式约束库容上下限、发电流量上下限 Aub, bub [], [] for t in range(T): row np.zeros(nvars); row[T t] 1 Aub.append(row); bub.append(V_max) row np.zeros(nvars); row[T t] -1 Aub.append(row); bub.append(-V_min) row np.zeros(nvars); row[t] 1 Aub.append(row); bub.append(Q_max) row np.zeros(nvars); row[t] -1 Aub.append(row); bub.append(-Q_min) # 等式约束水量平衡 V(t1) V(t) (Q_in - Q_out) * dt Aeq, beq [], [] for t in range(T): row np.zeros(nvars) row[T t] -1 # -V(t) row[t] dt # Q_out(t)*dt移到等式左侧 if t T - 1: row[T t 1] 1 # V(t1) beq.append(flow_in[t] * dt) else: # 末时刻强制 V(24) 回到初始库容体现日循环调度 beq.append(V_0 - flow_in[t] * dt) Aeq.append(row) res linprog(c, A_ubAub, b_ubbub, A_eqAeq, b_eqbeq, bounds[(0, None)] * nvars, methodhighs) print(求解状态:, res.message) Q_out res.x[0:T] V_opt res.x[T:2*T] # 风光出力风速/辐照转功率这里用模拟时序 wind_ms np.array([6,5,4,4,5,7,9,11,12,13,12,11, 10,8,7,6,6,5,4,4,3,3,2,2]) ghi np.array([0,0,0,0,0,50,200,450,700,850,900,880, 750,550,350,150,50,0,0,0,0,0,0,0]) p_wind_rated 50.0 # 风电装机 MW p_pv_rated 30.0 # 光伏装机 MW wind_cf 0.30 # 风电综合出力折算系数 pv_cf 0.80 # 光伏综合效率(温度/灰尘/逆变器) P_wind np.minimum(wind_ms / 12.0, 1.0) * p_wind_rated * wind_cf P_pv ghi / 1000.0 * p_pv_rated * pv_cf P_hydro coef_h * Q_out df pd.DataFrame({ Q_out: Q_out, V_mcm: V_opt, P_hydro_MW: P_hydro, P_wind_MW: P_wind, P_pv_MW: P_pv, total_MW: P_hydro P_wind P_pv }) print(df.head(10)) print(水电最大出力MW:, round(P_hydro.max(), 2))这段代码的逻辑说明决策变量分两段0到23是逐小时发电流量24到47是逐小时库容。水量平衡写成等式约束库容上下限和流量上下限写成不等式。把目标函数设为负的水电出力交给linprog做最小化风光出力不属于决策变量直接按折算系数算出来叠加。dt 3600 / 1e6是关键它把m3/s的流量换算成“每小时的百万m3”保证水量平衡两侧单位一致。末时刻等式把第24小时末库容拉回初始值80等价于“日循环调度”来水入库多少就发电用掉多少不留跨日余量。参数说明coef_h用密度1000、重力9.81、水头50米、效率0.85计算结果是0.417 MW/(m3/s)也就是每秒1立方米流量发0.417兆瓦。这里的水电最大出力约25MW和第3章案例里的50MW不需要完全一致实际项目里调度模型和财务模型要用同一套装机参数但教学演示可以分离。Q_max设到60m3/s是受水轮机过流能力限制不是越大越好如果把上限放太松汛期调度结果会偏乐观。注意这个最小模型刻意省略了弃水变量和坝前水位-出力耦合。日调节、来水小于装机过流能力的场景下够用一旦遇到汛期必须按第5.3节加弃水松弛变量否则库容约束直接冲突。3. 净现值计算把未来二十多年的发电收益折成今天的一笔钱3.1 净现值公式与现金流拆解哪些项必须进模型净现值NPV的标准公式是NPV -C_0 Σ CF_t / (1r)^tC_0是初始投资CF_t是第t年净现金流r是折现率N是运营期年限。水风光互补项目的现金流入几乎只有一项就是上网电费收入。现金流出包括运维成本、水资源费、保险费、大修基金、贷款还本付息。但可研阶段的NPV估算我一般只保留五个科目一是初始投资C_0包含机电设备、土建、接入系统、前期费用、建设期利息。水电站的单位千瓦造价明显高于风光但水电寿命长、出力稳定。二是年发电收入年发电量乘综合上网电价年发电量必须由第2章调度模型输出不要用装机容量乘年利用小时数粗算。三是年运维成本通常按初始投资的1%到2%估算风电场略高、光伏略低、水电介于中间混合项目要分类加总。四是发电衰减率光伏组件每年衰减约0.5%到0.8%风机有可利用率约束约97%水电站面临淤积导致出力递减NPV里不能不衰减。五是残值运营期末的设备残值与拆除清理费用。折现率取多少是一个决策问题而不是纯技术问题。国内可研阶段水电项目常用6%到8%风光项目因为风险略高常用7%到9%联合项目我一般取8%做基准然后做6%到10%的敏感性条带分析。折现率每抬1个百分点20年期的NPV大概缩水10%到15%这是项目决策里最敏感的单一参数。3.2 光伏风电水电三站联合NPV计算代码与基础参数调节import numpy as np # 基础参数 life 25 # 运营期25年 r 0.08 # 基准折现率8% tariff 0.38 # 综合上网电价元/kWh含绿证/现货溢价可调 # 初始投资单位亿元 capex { hydro: 4.5, # 50MW水电单位造价约9000元/kW wind: 2.8, # 50MW风电单位造价约5600元/kW pv: 1.5 # 30MW光伏单位造价约5000元/kW } C0 sum(capex.values()) # 年发电量拆分kWh/年由调度模型输出后替换这两个值 E_split {hydro: 1.2e8, wind: 0.8e8, pv: 0.6e8} # 运维成本占初始投资比例 opex_ratio {hydro: 0.015, wind: 0.020, pv: 0.010} opex_base {k: capex[k] * v for k, v in opex_ratio.items()} # 年衰减率水电最低光伏最高 degradation {hydro: 0.002, wind: 0.004, pv: 0.006} # 逐年现金流与NPV npv -C0 cashflows [] for t in range(1, life 1): f_h (1 - degradation[hydro]) ** (t - 1) f_w (1 - degradation[wind]) ** (t - 1) f_p (1 - degradation[pv]) ** (t - 1) E_t (E_split[hydro] * f_h E_split[wind] * f_w E_split[pv] * f_p) revenue E_t * tariff opex_t sum(opex_base.values()) * (1.005 ** (t - 1)) # 运维随通胀微涨 insurance revenue * 0.01 cf_t revenue - opex_t - insurance cashflows.append(cf_t) npv cf_t / (1 r) ** t print(f初始投资C0: {C0:.2f} 亿元) print(f第1年净现金流: {cashflows[0]:.2f} 亿元) print(f第25年净现金流: {cashflows[-1]:.2f} 亿元) print(fNPV(8%): {npv:.2f} 亿元) # 折现率敏感性6%~10% for rr in [0.06, 0.07, 0.08, 0.09, 0.10]: npv_rr -C0 for t in range(1, life 1): f_h (1 - degradation[hydro]) ** (t - 1) f_w (1 - degradation[wind]) ** (t - 1) f_p (1 - degradation[pv]) ** (t - 1) E_t E_split[hydro] * f_h E_split[wind] * f_w E_split[pv] * f_p revenue E_t * tariff opex_t sum(opex_base.values()) * (1.005 ** (t - 1)) cf_t revenue - opex_t - revenue * 0.01 npv_rr cf_t / (1 rr) ** t print(f折现率{rr:.0%}: NPV {npv_rr:.2f} 亿元)这段代码的逻辑说明把25年逐年现金流折现年电量考虑了三种电源的衰减率差异水电最低、光伏最高。运维成本引入1.005的年膨胀系数模拟人工和设备价格上涨。折现率敏感性循环直接输出六档NPV这个输出可以画成净现值-折现率曲线放进可研报告。基准条件下第1年净现金流大约0.84亿元第25年降到约0.73亿元NPV折现后约0.15亿元属于“刚过线”的项目。敏感性格局很有代表性折现率降到6%时NPV超过2.5亿升到10%就转负。参数说明tariff是最需要谨慎的参数。目前国内风光已进入平价阶段实际成交价在0.25到0.40元/kWh区间波动如果项目里有补贴、绿证交易或电力现货市场溢价要按“综合电价”计算而不是单看标杆电价。life设为25年是风光项目的常见运营期水电可以放到30年甚至40年但年限拉长后折现因子已经把远端现金流压得很薄30年以后的贡献对NPV影响很小。如果算出来NPV为负先别怀疑代码先检查电价假设和单位千瓦投资是否脱离当地实际。4. 不同网格装机容量与出力系数的敏感性分析网格尺度如何左右最优方案4.1 网格化资源评估为什么0.1°和0.5°的出力系数能差出10%“不同网格装机容量和出力系数”这个点核心是空间分辨率的敏感性。水风光资源评估通常使用再分析数据或卫星反演数据常见空间分辨率有0.05°约5.5km、0.1°约11km、0.25°约28km、0.5°约55km。分辨率越粗一个网格覆盖的地理范围越大网格内的风速、辐照、降水的空间异质性被平均掉。这个“平均掉”的影响是系统性的。风速分布是非线性的风机出力在切入风速到额定风速之间近似是风速的立方关系所以粗网格几乎必然高估风电出力系数。光伏辐照的均质化影响小一些但局部地形遮挡、山地云雾在粗网格里可能导致辐照偏差。水电的径流数据更依赖流域尺度网格大小关系到蒸散发和降水产流计算的汇流路径粗网格往往把洪峰削平汛期来水过程失真。因此做敏感性分析的标准动作是把同一套装机容量方案放到0.05°、0.1°、0.25°、0.5°四套网格上分别计算出力系数和年发电量统计差异。如果差异小于2%说明该区域资源空间均匀性好方案对这个参数不敏感如果差异超过10%常见于山地、海峡、复杂地形区域那么网格分辨率应该被当作项目可行性的不确定因子在NPV里做概率化处理而不是用单值。4.2 网格扫描与装机容量配置联动代码实现import numpy as np import pandas as pd # 模拟四套网格下的资源评估输出 # 实际项目中这里读的是不同分辨率的nc/tif文件这里用模拟数据演示扫描逻辑 grids [0.05deg, 0.10deg, 0.25deg, 0.50deg] res_assess { 0.05deg: {wind: 7.8, ghi: 1650, flow: 38.5}, 0.10deg: {wind: 8.0, ghi: 1640, flow: 38.8}, 0.25deg: {wind: 8.4, ghi: 1610, flow: 39.6}, 0.50deg: {wind: 8.9, ghi: 1570, flow: 40.8} } # 候选装机方案: (水电MW, 风电MW, 光伏MW) capacity_plans [ (40, 45, 25), (50, 50, 30), (50, 60, 40), (60, 70, 50), ] # 出力系数换算可研近似版 def capacity_factor(wind_ms, ghi_kwh, flow, plan): hydro_mw, wind_mw, pv_mw plan # 水电流量相对设计流量46m3/s cf_h min(flow / 46.0, 1.0) * 0.70 # 风电简化功率曲线4~12m/s线性区 cf_w max(0, min((wind_ms - 4) / (12 - 4), 1.0)) * 0.95 # 光伏1kWp年发电量≈GHI*PRPR0.80 cf_p ghi_kwh * 0.80 / 8760.0 return cf_h, cf_w, cf_p results [] for grid in grids: data res_assess[grid] for plan in capacity_plans: cf_h, cf_w, cf_p capacity_factor( data[wind], data[ghi], data[flow], plan) hydro_mw, wind_mw, pv_mw plan e_h hydro_mw * cf_h * 8760 / 1e3 # GWh e_w wind_mw * cf_w * 8760 / 1e3 e_p pv_mw * cf_p * 8760 / 1e3 total_e e_h e_w e_p results.append({ grid: grid, plan: str(plan), cf_h: round(cf_h, 3), cf_w: round(cf_w, 3), cf_p: round(cf_p, 3), GWh: round(total_e, 1) }) df pd.DataFrame(results) pivot df.pivot_table(indexplan, columnsgrid, valuesGWh) print(pivot) # 每个方案在不同网格下的最大相对偏差 def max_grid_dev(s): return (s.max() - s.min()) / s.mean() dev df.groupby(plan)[GWh].agg(max_grid_dev) print(dev.round(3))这段代码的逻辑说明先模拟四套分辨率的资源评估结果再对四种装机容量方案逐一计算水电、风电、光伏出力系数和年发电量用pivot输出“方案×网格”矩阵最后统计每个方案在不同网格尺度下年发电量的最大相对偏差。模拟数据里隐含了一个典型趋势网格越粗风速越高、辐照略低、径流略大。这个趋势在某些区域会反转比如孤立山峰区域粗网格反而低估风速所以敏感性扫描不能省。参数说明capacity_factor函数里的系数是示意性的实际项目中应该替换为风电机组厂商提供的功率曲线插值、光伏组件实测PR值和电站设计水头下的水轮机效率曲线。光伏的0.80是PR系统综合效率覆盖温度、灰尘、逆变器损失年辐照1650kWh/m2的地区1kWp组件年发电约1320kWh对应出力系数0.151这是中国西北地区的常见数值。严格的风电出力系数应该用小时间隔的风速序列逐点代入功率曲线再求平均用年均风速代入功率曲线只是可研阶段的近似粗网格下这个近似会额外放大偏差。如果max_grid_dev超过5%这个方案对网格分辨率敏感在可研阶段需要加密网格或改用更精细的局部资源数据。5. 水风光互补调度的常见问题与避坑记录从数据到求解器的五个坑5.1 径流、风速、辐照三条时序的时间基准错位现象调度模型跑出来的出力曲线在小时尺度上断崖式跳变电量加总比手工估算低15%以上。原因三个数据源的时间戳不是同一个时区。风速用了UTC时间辐照用的是地方太阳时径流则是日平均数据插值成小时。三者错位3到5个小时互补调度的“补位”逻辑完全失效水电在风光低谷时段的库容分配全错。解决在项目启动第一天做统一时间基准。先把三条数据重采样到同一网格和同一时区再检查跨日、跨月、夏令时边界。我习惯用数据源自带的原始时间戳先对齐到UTC再统一转换为地方时。时间基准这类错误很难靠肉眼发现有一条简单检查方法把水风光三条出力曲线叠加画出来如果水电负荷变化和风光低谷曲线错位超过2小时优先怀疑时间基准。5.2 出力系数口径不统一容量基准一个天上一个地下现象同一个风电场两份方案书里写的出力系数分别是0.28和0.36但上网电量完全相同。原因一个用额定容量做分母一个用“铭牌容量×可利用率”做分母另一个把限电电量也算进分子。口径不一致导致NPV和调度模型两个环节各自使用不同系数最终经济评价出现矛盾。解决在代码里固定出力系数定义CF 实际发电量 /额定装机 × 8760小时。所有网格、所有方案的出力系数一律用这个口径计算限电和检修停机单独记录不进分子。这个约定要在项目一开始就写进数据处理脚本里而不是等到出报告时再统一。5.3 弃水约束做错库容直接顶到上限调度解不靠谱现象汛期时段水电出力接近满发但库容越界线性规划报“无可行解”或者水库水位恒定在最高位完全不落调度曲线看起来异常平滑。原因水量平衡方程里漏了弃水变量导致库容上限约束在来水大于发电用水时必然冲突。因为线性规划里不能写“如果库容满则弃水”只能用变量表达很多初学者把Q_out限制为Q_in和V_max的函数这在优化模型里根本表达不了。解决增加弃水松弛变量Q_spill并放进水量平衡等式V(t1) V(t) (Q_in - Q_out - Q_spill) × dt给弃水变量加一个极小的成本系数让求解器在“弃水”和“越库”之间自动选择。伪代码如下# 在2.2节的变量段上扩展第三段: [Q_out_t, V_t, Q_spill_t] # 水量平衡变为: # V(t1) V(t) (Q_in(t) - Q_out(t) - Q_spill(t)) * dt # 目标函数中给Q_spill加大惩罚系数例如在c向量末尾追加 0.01加了这个变量之后汛期库容不会再顶到上限弃水会集中在来水峰值时段这符合水电站实际调度规律。5.4 折现率装的是“拍脑袋值”敏感性分析只做了一档现象NPV计算结果在评审会上被质疑折现率从8%调到9%项目就从正收益变成负收益决策悬空。原因折现率取了单一固定值没展示6%到10%区间的NPV变化。水风光联合项目叠加了电改政策、市场电价、来水偏枯三重不确定性单点NPV不具备决策参考价值。解决NPV输出处顺手生成折现率敏感性表同时做概率化分析。抽取来水、风速、辐照、电价四个随机变量跑几百次蒙特卡洛画出NPV直方图比单点计算有说服力得多。敏感性代码不用额外库第3章里那个for循环已经够用。5.5 网格分辨率混用资源数据0.25°、装机数据在乡镇边界上硬算现象网格敏感性分析结果异常相邻两个方案的年发电量差异高达20%明显不符合物理规律。原因一套项目里同时用了0.25°的风速数据、0.1°的辐照数据和流域水文站点的径流数据三条数据叠加之后网格边界根本无法对齐空间插值把伪信号当成了真差异。解决先做资源数据的网格归一到统一分辨率。推荐用面积加权聚合把细网格归到粗网格或用双线性插值把粗网格加密。复杂山地用双线性不如最近邻稳健因为线性插值会把山顶和峡谷的值平滑掉。做完归一化再做敏感性扫描同时把每个网格内站点数量和插值方法记录在元数据里方便审计。6. 进阶验证用LCOE和内部收益率交叉确认NPV结论NPV算出来为正不代表这个方案在同等风险下足够有竞争力。评审时通常会追问两个衍生指标内部收益率IRR和平准化度电成本LCOE。IRR就是让NPV等于0的折现率LCOE是全生命周期成本除以全生命周期发电量单位元/kWh。NPV告诉绝对收益LCOE告诉成本竞争力IRR告诉风险的耐受边界。建议在调度-经济一体化模型跑通后把下面这段验证代码加进去。from scipy.optimize import brentq def calc_lcoe(c0, opex_t_list, e_t_list, r): 平准化度电成本折现成本 / 折现电量 cost_pv c0 sum(opex / (1 r) ** t for t, opex in enumerate(opex_t_list, start1)) e_pv sum(e / (1 r) ** t for t, e in enumerate(e_t_list, start1)) return cost_pv / e_pv def calc_irr(c0, cf_list): 内部收益率brentq在(0.1%, 100%)找NPV0的根 f lambda r: -c0 sum(cf / (1 r) ** t for t, cf in enumerate(cf_list, start1)) try: return brentq(f, 0.001, 1.0) except ValueError: return np.nanLCOE算出来后和当地标杆电价或市场均价比较如果LCOE高于电价1.5倍以上说明发电成本过高即使NPV因为某些补贴为正也属于靠政策养着的项目。IRR扫描时如果发现8%折现率下NPV为正、9%就转负说明项目风险容忍度很低要主动降低装机容量或延长运营期而不是调高出力系数去迎合评审。我的习惯是给每个方案建立一张“三指标卡”NPV、IRR、LCOE并列同时标注网格敏感性的最大偏差值。如果最大偏差超过5%三指标卡上加一行“资源数据分辨率需加密”的批注。这套检查做完前面那些模棱两可的参数才真正立得住。希望帮到你。本文还有配套的精品资源点击获取