简介基于鲸鱼优化算法WOA与CVaR风险度量的售电公司购售电优化模型MATLAB实现面向电力市场交易决策人员、电力系统/优化算法方向研究生及工程师解决售电公司在多重购电渠道与两类售电合同下的最优购电策略与收益风险平衡问题。模型完整覆盖购电侧长期市场、现货、可再生能源、分布式电源及储能租赁五类业务售电侧含均一电价与实时电价合同并给出基于CVaR的风险评价方法。压缩包共19个文件以13个m脚本和4个mat数据文件为主体另含1个txt说明与1份12000字论文doc总计约501KB代码与论文对应清晰便于直接运行复现。已有568人学习下载通过WOA的包围、螺旋捕食及随机觅食三个步骤对购售电收益风险目标函数寻优最终获得并验证最优购电策略适合作为交易策略优化、风险建模及智能算法应用的完整参考。1. WOA与CVAR在售电公司购售电策略中的定位售电公司做购电侧和售电侧模型真正要解决的其实是一个组合优化问题在批发市场买什么、买多少在零售市场以什么价卖给用户、签哪种合约同时要控制极端行情下的亏损风险。传统的期望值优化只盯着平均收益但电力现货价格分布是厚尾的对价格尖峰无能为力。CVAR条件风险价值正好补上这个缺口——它度量的是“最差那5%或1%场景下的平均损失”在收益和风险之间给出显式平衡。而带CVAR的购售电模型目标函数通常非凸、多峰梯度下降类求解器容易陷在局部解里这就是引入WOA鲸鱼优化算法的动机WOA是一种无梯度群体智能算法不需要目标函数连续可导对约束边界、非线性惩罚项甚至不可微的0-1变量都有不错的适应度。本文直接沿着“CVAR建模→WOA求解→MATLAB实现→结果验证”这条线展开适合正在做电力市场方向本科或硕士论文的学生以及想在购售电决策里引入风险控制量化手段的工程师。2. 购电侧CVAR模型构建风险度量方式与场景生成2.1 CVAR的数学定义与在购电组合中的含义CVAR的定义依赖于VAR风险价值。设购电组合某一场景下的成本为随机变量 (f(x, \xi))其中 (x) 是决策变量购电量、购电结构(\xi) 是随机场景现货价格、新能源出力、负荷水平那么在置信水平 (\alpha)常用0.95下VAR是满足 (P(f(x,\xi) \text{VaR}_\alpha) \le 1-\alpha) 的分位数CVAR则是超过该分位数的条件期望[ \text{CVaR}\alpha(x) E\left[f(x,\xi) \mid f(x,\xi) \text{VaR}\alpha(x)\right] ]在购电侧模型里f代表总购电成本。与VAR相比CVAR具备次可加性和凸性它不会因为把组合拆开就“变安全”所以用CVAR作为风险约束比用VAR更稳健。实际建模中是把它写成线性规划可解的形式借助场景集近似[ \text{CVaR}\alpha(x) \approx \text{VaR}\alpha \frac{1}{(1-\alpha) S} \sum_{s1}^{S} \max\left( f(x, \xi_s) - \text{VaR}_\alpha, 0 \right) ]这里 (S) 是场景总数。MATLAB实现时不需要显式求解这个分位数而是引入辅助变量 (z_s) 来线性化 (\max) 项。2.2 购电侧模型包含的三个典型模块常见做法是让购电侧由三部分构成中长期双边合约电量、日前现货市场购电、新能源风/光出力。三者的成本特性完全不同。中长期合约价格锁定但灵活性差现货市场价格波动剧烈但可以按需调整新能源边际成本近乎零但出力不确定性强需要风控来保证收益稳定。购电成本表达式如下[ C_{\text{total}} \sum_{t1}^{T}\left( C_{\text{contract},t} \cdot Q_{\text{contract},t} \lambda_{\text{spot},t} \cdot Q_{\text{spot},t} C_{\text{renew},t} \cdot Q_{\text{renew},t} \right) ]其中 (Q_{\text{spot},t}) 是需要实时决策的变量受出力上限约束(Q_{\text{contract},t}) 通常是固定值或按周调整。MVAR模型在这个基础上加入一个风险约束[ \text{CVaR}\alpha\left(C{\text{total}} \beta \cdot Q_{\text{punish}} - R_{\text{expected}}\right) \le \rho ](\rho) 是风险容忍上限(\beta \cdot Q_{\text{punish}}) 是偏差惩罚项用来限制实际成本偏离预期成本的幅度。2.3 场景生成用历史现货价格做蒙特卡洛模拟要算CVAR必须有场景集。一般先用历史现货价格拟合分布常见做法是用几何布朗运动GBM模拟下一周期的价格路径或者直接拿去年同期数据做keyboard bootstrap。MATLAB里GBM模拟代码如下% 场景生成几何布朗运动模拟现货价格 % P0: 初始价格mu: 漂移率sigma: 波动率dt: 时间步长nSteps: 步数nScen: 场景数 rng(42) P0 350; % 初始现货价格单位: 元/MWh mu 0.002; % 日漂移率由历史数据统计得到 sigma 0.35; % 日波动率同样由历史数据估计 dt 1; % 单位时间设为1天 nSteps 30; % 30天购电周期 nScen 500; % 500个随机场景 % 预分配场景矩阵每一行是一条价格路径 Sce zeros(nScen, nSteps1); Sce(:,1) P0; for s 1:nScen for t 1:nSteps Sce(s, t1) Sce(s, t) * exp((mu - 0.5*sigma^2)*dt sigma*sqrt(dt)*randn()); end end % 截掉起始价格列只保留未来30天 Sce Sce(:, 2:end);说明一下参数含义漂移率mu决定价格长期趋势波动率sigma决定价格路径的离散程度这两个值必须用同一时间窗口的历史数据来标定不能用拍脑袋值。nScen500是经验值场景太少会让尾部概率估计不稳太多则计算代价上升。生成场景后建议先画一下场景分布直方图验证模拟价格不会出现大量负值——GBM对极端波动加负漂移时可能出现非物理负价负价场景需统一截断为0。2.3.1 风电出力的场景耦合问题购电侧模型中的新能源出力与现货价格之间存在相关性不能独立建模。较稳妥的做法是直接对“风功率-现货价格”做二维Copula采样而不是分别做一维模拟再乱序拼接。如果手里只有历史数据且嫌Copula实现繁琐可以退一步使用正向抽样法先抽样风功率场景再在当前风功率条件下近似估算价格的条件分布。代价是模型的边缘分布会略失真但对于以CVAR为核心的论文和项目验证来说足够用。2.4 购电成本的向量化计算把场景矩阵直接与决策变量做矩阵乘法避免用for循环累加成本这会显著提升MATLAB运行速度% 假设q_contract和q_renew是已知向量q_spot是待决策变量 % Sce是nScen x T的价格场景矩阵 % 汇总每个场景下的总购电成本向量 (nScen x 1) C_contract q_contract * c_contract; % 中长期合约总成本标量 C_spot sum(Sce .* repmat(q_spot, nScen, 1), 2); % 每个场景的现货购电成本 C_renew q_renew * c_renew; % 新能源购电成本标量因为价格固定 C_total_vec C_spot C_contract C_renew; % 每个场景的总成本repmat将决策变量横向扩展到与场景矩阵相同维数然后在第二维上求和。这样写不仅代码短而且能利用MATLAB内置的多线程计算。对于500个场景、30个时间阶段这个矩阵乘法在普通笔记本上也是毫秒级。3. 鲸鱼优化算法WOA寻优核心循环与适应度函数设计3.1 WOA算法为何能匹配CVAR这类非光滑目标CVAR目标函数通过线性化辅助变量后虽然理论上可导但在MATLAB实际实现中往往还带有0-1变量比如是否启用某个合约、阶梯价格区间、风电偏差惩罚项这些让目标函数布满平台区和跳变点。梯度类方法需要对梯度做次梯度近似调起来非常痛苦。WOA的搜索方式是位置向量直接更新完全不需要计算梯度天然适合这类黑盒目标函数。WOA的机理如下每头鲸鱼代表候选解位置向量对应一组购电量决策。算法分三个环节——包围猎食、气泡网攻击、随机搜索。包围猎食阶段是当前最优解吸引气泡网攻击阶段用对数螺旋收缩来精细搜索随机搜索则保持多样性、跳出局部最优。3.2 适应度函数如何嵌入CVAR适应度函数有两个任务既要让总购电成本尽量低又要控制尾部损失。将约束条件用惩罚函数形式并入目标代码如下function fitness woa_objective(q_spot, Sce, params) % 输入 % q_spot : 长度为T的现货购电量决策向量 % Sce : nScen x T 的现货价格场景矩阵 % params: 包含合约电价、风险系数beta、置信度alpha等 % % 输出 % fitness: 标量包含成本期望、CVAR和惩罚项的综合指标 T params.T; nScen size(Sce, 1); % 1) 计算每个场景的总购电成本 C_contract params.contract_price * params.q_contract; % 固定 C_renew params.renew_price * params.q_renew; % 固定 C_spot_vec (Sce .* repmat(q_spot(:), nScen, 1)) * ones(T, 1); cost_vec C_spot_vec C_contract C_renew; % 2) 按升序排列计算VaR和CVaR sorted_cost sort(cost_vec); alpha_idx max(1, floor(params.alpha * nScen)); % 0.95分位索引 VaR_val sorted_cost(alpha_idx); tail_mean mean(sorted_cost(alpha_idx:end)); % 尾部平均 tail_loss max(0, tail_mean - VaR_val); % 相对风险值 % 3) 期望成本 CVAR风险项 爬坡约束惩罚 expected_cost mean(cost_vec); fitness expected_cost params.risk_coef * tail_loss; % 4) 爬坡约束惩罚现货购电量不允许跳变过大 ramp_penalty sum(max(0, abs(diff(q_spot)) - params.ramp_limit)); fitness fitness params.ramp_coef * ramp_penalty;这段代码里有两个容易踩的细节。一是CVAR计算用排序法而不是线性规划形式原因是在适应度函数里做线性化辅助变量会引入大量额外维度让WOA的搜索空间维度从T变成TnScen收敛速度明显变慢排序法对500个场景完全可接受且结果一致。二是爬坡约束用二次惩罚而非硬约束把违规量加进目标函数这样WOA迭代初期会优先修正严重越限的解。3.3 WOA主循环的MATLAB实现WOA自身的实现比较标准主循环包括种群初始化、迭代更新三套位置更新公式。下面是核心循环代码function [best_sol, best_fit, convergence] woa_solver(obj_func, dim, lb, ub, nWhales, maxIter) % obj_func : 适应度函数句柄 % dim : 决策变量维度 (T) % lb, ub : 决策变量下上和上界向量 % nWhales : 种群大小 % maxIter : 最大迭代次数 % 初始化位置LHS拉丁超立方采样比随机均匀分布更均匀 X lhsdesign(nWhales, dim) .* (ub - lb) lb; fitness zeros(nWhales, 1); for i 1:nWhales fitness(i) obj_func(X(i,:)); end [best_fit, best_idx] min(fitness); best_sol X(best_idx, :); convergence zeros(maxIter, 1); for iter 1:maxIter a 2 - 2 * iter / maxIter; % 线性衰减参数控制收敛 for i 1:nWhales r1 rand(); r2 rand(); A 2 * a * r1 - a; C 2 * r2; p rand(); b 1; % 螺旋形状常数论文里通常取1 l (rand() - 0.5) * 2; % [-1, 1]之间 if p 0.5 if abs(A) 1 % 包围捕食向当前最优鲸鱼移动 D abs(C .* best_sol - X(i,:)); new_pos best_sol - A .* D; else % 随机搜索随机选一条鲸鱼摆脱局部极值 rand_idx randi(nWhales); X_rand X(rand_idx, :); D abs(C .* X_rand - X(i,:)); new_pos X_rand - A .* D; end else % 螺旋气泡网更新对数收缩 D abs(best_sol - X(i,:)); new_pos D .* exp(b .* l) .* cos(2*pi*l) best_sol; end % 边界处理 更新适应度 new_pos max(min(new_pos, ub), lb); new_fit obj_func(new_pos); if new_fit fitness(i) X(i,:) new_pos; fitness(i) new_fit; end end [current_best, idx_best] min(fitness); if current_best best_fit best_fit current_best; best_sol X(idx_best, :); end convergence(iter) best_fit; end end参数含义要看清nWhales是鲸鱼数量大一点搜索广但耗时长a从2线性降到0是WOA的收敛机制前期全局搜索后期局部精细搜索。p 0.5且|A| 1时执行包围捕食|A| 1时强制随机搜索这能有效防止种群过早集中到某个不合理解。这里每一维都独立计算随机数实际效果比共享同一个随机数好。3.4 与粒子群PSO的性能对比WOA的适用边界WOA不是万能的。当决策变量维度低但约束高度严格时PSO带惯性权重通常收敛更快当目标函数存在大量局部极值且维度较高时WOA的螺旋搜索和随机搜索组合更有优势。对购电侧模型决策变量维度就是时间阶段数TT在24到168之间维度中等偏高。可以做一个简单的基准测试脚本对比二者但重点观察收敛点而不是收敛速度——找到的全局最优解质量比前几十次迭代的下降速度重要得多。另外必须注意的是WOA的标准形式连续分布算法对整数变量如合约份数需要做取整处理或者用TOPK的映射方式。购电侧模型中如果出现“合约数量是整数”这种约束建议在目标函数内部取整并加微小惩罚不要直接对位置向量取整后在主循环中反馈那样会让边界处的搜索失效。4. 售电侧模型与购电侧闭环联动定价策略和用户分类4.1 售电侧模型的分时电价与三类用户购电侧确定下来以后售电侧的任务是制定零售价格套餐让售价覆盖购电成本且留有利润。现实中售电公司一般把用户分成三类。工业用户峰谷差明显但对电价敏感商业用户负荷平稳且价格接受度高居民用户数量多、单体负荷小。分时电价策略为 ( P_{\text{retail},t} P_{\text{base}} \Delta P_t )其中 (\Delta P_t) 是时段调整量分峰、平、谷三档。售电侧的优化目标是在给定用户需求价格弹性的条件下最大化售电利润并控制售电侧风险。利润函数考虑动态需求响应[ R_{\text{sale}} \sum_{t1}^{T}\left( P_{\text{retail},t} \cdot Q_{\text{demand},t}(P_{\text{retail},t}) \right) ]其中 (Q_{\text{demand},t}) 是用户在当前电价下的用电量用线性弹性模型近似% 用户负荷需求与电价的线性弹性模型 % 基础负荷 Q0弹性系数 e一般为负基准电价 P0调整后电价 P function Q_demand demand_response(P_retail_t, Q0, e, P0) % 弹性系数含义: e-0.15 表示价格上升10%需求下降1.5% Q_demand Q0 .* (1 e .* (P_retail_t - P0) ./ P0); end4.2 从购电成本到售电价差闭环目标函数把购电侧模型和售电侧模型联通的关键是价差。每个场景下的购电成本是 (C_{\text{spot}}(\xi_s))而售电收入由定价策略决定两者之间存在相关性——现货价格高时如果售电价没跟上利润就被压缩。闭环目标函数应直接定义净利润场景向量function overall_fitness integrated_obj(price_retail_vec, q_spot_vec, Sce_p, Sce_load, params) % 输入两个决策数组售电套餐价和现货购电量 % Sce_load 是各场景下的负荷需求矩阵 nScen size(Sce_p, 1); profit_vec zeros(nScen, 1); for s 1:nScen Q_load demand_response(price_retail_vec, params.Q0, params.elasticity, params.P0); revenue sum(price_retail_vec .* Q_load); cost sum(Sce_p(s,:) .* q_spot_vec) params.fixed_cost; profit_vec(s) revenue - cost; end % 净利润用期望 CVaR风险惩罚 exp_profit mean(profit_vec); sorted_profit sort(profit_vec); alpha_idx max(1, floor((1 - params.alpha) * nScen)); VaR_loss sorted_profit(alpha_idx); cvar_loss mean(sorted_profit(1:alpha_idx)); % 左尾平均即最差收益的平均 fitness -(exp_profit - params.risk_coef * abs(cvar_loss)); end注意这里风险度量的方向变了购电侧最小化成本时用CVAR约束上尾售电侧最大化利润时CVAR度量的是利润分布左侧尾部——也就是最差那一批场景的平均利润。所以代码里排序后取的是前alpha_idx个而不是后alpha_idx个。方向搞反会让模型变成“鼓励极端亏损”。4.3 WOA在两个模块之间如何衔接最终模型中WOA的决策向量应当是[售电价格方案, 现货购电量]的拼接向量而不是分开跑两个WOA。分开优化必然导致次优售电侧最优定价是按某个固定购电成本做的而购电侧优化又不知道售价的约束。将两组决策拼接在同一位置向量中让WOA同时搜索才是二个模块联合优化的正确做法。决策向量长度变为 (2T)维度增加后需要适当增大种群数量。经验值表格如下决策维度种群数最大迭代数CVAR置信度风险系数区间T24301000.950.1 ~ 0.5T48401500.950.1 ~ 0.5T168603000.9 或更低0.05 ~ 0.3风险系数取多少要实际跑数据看。系数太大模型会明显降低期望利润来换取波动率下降系数太小CVAR的影响几乎可以忽略模型等价于确定性优化。论文写作时普遍的做法是给出风险系数从0到1的横扫曲线说明不同风险偏好下的“利润-CVAR前沿”而不是只给一个点。5. 结果验证与调参技巧从代码到论文结论的衔接5.1 收敛曲线和场景稳定性检验WOA是随机算法必须检验多次运行的一致性。做法是固定随机种子跑100次记录每次的最优值和平均收敛曲线计算最优解的标准差。标准差与均值的比值如果超过10%说明种群数太小或迭代次数不足需要加大。收敛曲线不应该只看单次要看95%置信区间带。5.2 CVAR置信度与场景数的敏感性分析这是论文中很有分量的一个实验。分别取置信度0.80、0.85、0.90、0.95、0.99观察最优解的变化规律。通常置信度越高最优解对应的现货购电比例越低合约电量占比越高这是符合直觉的——高置信度要求更强对冲。场景数从100递增到2000验证CVAR值是否稳定如果2000个场景与500个场景得到的CVAR值差超过5%就说明场景生成过程有问题通常要么是价格模拟参数没标定好要么是风电场景没与价格场景正确耦合。5.3 论文结构组织和图表选型12000字论文的基本框架怎么与代码对应模型章节展示3.2节的目标函数和约束条件算法章节展示3.3节WOA伪代码和参数表算例章节应当放三个结果图——收敛曲线对比WOA vs PSO、利润-CVAR前沿图和购电结构分时柱状图。这三个图覆盖了“算法性能模型效果决策可视化”三个审稿人最关注的维度。注意图不要用MATLAB默认的jet配色改用parula或手动指定色系在黑白打印下也能分辨这在论文盲审阶段是加分项。5.4 收敛速度优化技巧大规模场景下WOA速度瓶颈在适应度函数而不是WOA本身。如果场景数达到5000以上应把适应度函数里耗时的场景循环改为完全向量化。一个常用技巧是预计算价格场景矩阵与决策无关的部分每次调用函数时只计算变化的部分。另一个实用技巧是在迭代后半段把场景数从全量减到千分之一子样本每40代做一次全量评估。这种“粗评估定期精评”策略在实测中能把总运行时间压缩到原来的三分之一而最优解几乎不受影响。5.5 如果结果不理想优先检查的方向运行模型后如果发现WOA找到的“最优解”违反爬坡约束或CVAR约束首选检查惩罚项的系数是否覆盖了目标值量级。惩罚系数设得太小时约束违规在城市比成本更便宜求解器就会选择违规路径。应该先把只含惩罚项的目标单独跑一遍确认默认解的违规量为零再叠加成本项。其次检查变量边界是否合理现货购电量上限通常取负荷峰值的30%到50%或者按照电网公司与售电公司签订的购电合同上限设置不能随意取一个很大的数否则搜索空间里有大量物理上不合理的区域算法浪费在无效区域里。本文还有配套的精品资源点击获取