简介面向2021年电工杯数学建模B题这份资料将赛题思路解析与案例程序打包成套适合备赛选手、指导教师以及所有希望提升数学建模与电气工程应用能力的学习者。资源共包含70个文件压缩后约22.71MB其中37个mat数据文件提供仿真实测数据12个m脚本承载模型求解与绘图代码另有tif/png可视化图片、xlsx中间表格、pdf/doc完整说明文档等便于按数据、算法、结果、报告四类快速定位。目前已有3907人学习浏览。内容上既有对B题背景的深入分析和模型构建思路也有可直接运行调试的Matlab/Python案例程序初学者可借助注释代码理解从实际问题到数学模型的转化进阶者则可参考其算法设计并在电力系统优化、故障诊断等场景中复用。通过研读这些程序与数据读者能够较全面地掌握一类工程建模题的完整解法为后续参赛或科研积累可迁移的实操经验。1. 一份电工杯B题案例程序专治赛期“优化不收敛”2021年电工杯数学建模B题里的微电网经济调度是那种拿到题能列式子、跑起来却反复翻车的典型优化问题。负荷、光伏、风电、分时电价四类数据一进来变量维数、单位换算、爬坡约束和储能荷电状态会同时较劲很多队伍卡在“花了半天把约束写全粒子群一跑却不收敛”。这套案例程序把整个任务拆成数据读取、目标函数、粒子群迭代、出图四个独立模块用MATLAB在几分钟内给出一组可行的24小时调度方案。下面从程序结构开始拆讲清楚怎么跑、怎么改。2. 程序框架与解题思路从赛题描述到48维决策向量2.1 赛题要做的事微电网日前经济调度题目背景就是一个装了柴油发电机、储能、光伏和风电的并网型微电网24小时负荷曲线已知光伏和风电出力按典型日给出分时购电价也定好了。要做的决定是每个小时柴油机发多少、储能充还是放、从电网买多少或卖多少最后让一整天的总运行成本最低。把问题写成数学形式是这样决策变量PG(t)、PB(t)t1,…,24 目标min Σ_t [ a·PG(t)² b·PG(t) c π_buy(t)·max(Pnet(t),0) - 0.8·π_buy(t)·max(-Pnet(t),0) ] 其中 Pnet(t) Pload(t) - Ppv(t) - Pwind(t) - PB(t) 约束PG_min ≤ PG(t) ≤ PG_max |PG(t)-PG(t-1)| ≤ ramp |PB(t)| ≤ ESS_Pmax SOC_min ≤ SOC(t) ≤ SOC_max决策向量拼接成48维前24维是柴油机每小时出力后24维是储能功率。这里有个贯穿全套代码的符号约定——PB大于0表示储能充电等于从负荷侧吸收功率PB小于0表示放电向负荷供电。买电和卖电的价格不等卖电按购电价的80%计算贴近分布式光伏“自发自用、余电上网”的场景。这个模型里有几个容易写错的地方。第一SOC递推不是线性公式充电效率和放电效率不一样放电时要做除法而不是乘法第二爬坡约束是相邻小时之间的差分第1小时没有历史值可以不约束或单独处理第三如果题目要求储能日内复用即每天结束SOC回到初始值那就得在约束或目标里补上“两端守恒”的惩罚项。从竞赛角度看这道题用线性规划也能解Pnet、PG、PB都是连续变量目标函数里带max的部分可以拆成购、售两个非负变量线性化。案例程序之所以选粒子群一方面考虑到不少参赛队对线性规划建模不熟另一方面PSO写起来直观、调参空间大评阅时也能清晰看到模型思考过程不至于被一套黑盒求解器盖住工作量。2.2 解压后的目录结构数据、代码、文档各管一摊整个压缩包解开之后是下面这个布局B_case/ ├── data/ │ ├── load24.csv │ ├── pv24.csv │ ├── wind24.csv │ ├── price24.csv │ └── ESS_param.xlsx ├── code/ │ ├── main.m │ ├── read_data.m │ ├── objective_fun.m │ ├── pso_solve.m │ └── plot_schedule.m ├── result/ │ └── 运行后生成的图和调度数据 └── doc/ └── B题题目说明.pdfdata目录放原始数据code目录放全部脚本result目录在程序运行后自动生成输出图。这类结构在竞赛资源包里很常见核心思想是让数据、逻辑和输出分离——比赛后期换一组数据跑批只改data目录里的CSV代码一行不用动。程序执行顺序是main.m → read_data.m → objective_fun.m → pso_solve.m → plot_schedule.m。read_data返回带全部参数的结构体Dobjective_fun被pso_solve作为函数句柄反复调用最后统一出图。注意所有code里的脚本都用../data/这种相对路径引用原始数据所以别把某个子目录单独拷到别处运行哪怕只是换了一级目录名都可能直接报“找不到文件”。从这个意义上讲解压后保持目录结构原样就是第一个运行前提。main.m的主流程值得先看整个调用关系都在里面% main.m: 一键运行微电网日前经济调度 clear; clc; close all; D read_data(); % 读取24小时基础数据 dim 24 * 2; % 决策变量维数24h柴油出力 24h储能功率 lb [ones(24,1)*D.para.PG_min; ones(24,1)*(-D.para.ESS_Pmax)]; % 下限储能最小功率为负 ub [ones(24,1)*D.para.PG_max; ones(24,1)*D.para.ESS_Pmax]; % 上限储能最大功率为正 opt struct(swarm_size, 100, ... % 种群规模 max_iter, 500, ... % 最大迭代次数 c1, 2.0, c2, 2.0, ... % 个体学习因子 / 全局学习因子 w0, 0.9, w1, 0.4); % 惯性权重的起止值 [x_best, f_best, curve] pso_solve((x) objective_fun(x, D), ... dim, lb, ub, opt); plot_schedule(x_best, D, schedule_result.png); disp([最优成本 , num2str(f_best)]);这个脚本把变量边界直接写在下界/上界向量里而不是塞进目标函数做硬编码这是调参时最省事的写法。后面想加第三类决策变量比如可中断负荷比例只需要把dim改成24*3再在lb和ub里补一段0到1的边界目标函数同步读取就行。2.3 技术栈取舍为什么这里用MATLAB而不是Python很多参赛队第一反应是用Python写赛后复盘也觉得Python库多、更通用但这个包的实现选的是MATLAB原因很实际。一是MATLAB的矩阵语法和这个题目完美匹配。PG、PB、Pload这些全是24维向量mean、max、diff这类运算直接作用在数组上目标函数十几行就能写完。二是绘图方便plot、stairs、area几行就能出一张符合论文格式的图比赛修改比Python里调matplotlib省心。三是评阅老师对MATLAB这套语法太熟悉读代码成本低。四是比赛期间网络不一定稳定Python环境哪怕只是装个pandas和matplotlib也可能碰上镜像源问题MATLAB解压即用省掉一个潜在变量。如果包里有一份Python脚本我一般建议放在“数据初检”这个角色上比如提交前用pandas快速扫一遍CSV里有没有空值和数量级异常再回MATLAB正式跑。正式求解和出图别跨语言切换不然调试时来回换环境非常痛苦。3. 核心代码拆解从CSV到调度曲线的一条链路3.1 数据读取与预处理read_data.mread_data负责把五个文件读进来顺便做一层轻量清洗。CSV原始数据都是24行一列纯数值无表头readmatrix直接读就能拿到24x1向量。若文件带表头或单位注释行readmatrix在R2019a之后会自动把表头识别成字符串行跳过这点在换自己数据时要注意。function D read_data() % read_data 读取微电网24小时数据返回一个结构体D % D.load / D.pv / D.wind / D.price 均为24x1列向量 % D.para 保存柴油机与储能参数 base ../data/; D.load readmatrix([base, load24.csv]); % 24x1负荷功率 MW D.pv readmatrix([base, pv24.csv]); % 24x1光伏出力 MW D.wind readmatrix([base, wind24.csv]); % 24x1风电出力 MW D.price readmatrix([base, price24.csv]); % 24x1购电价 元/MWh % 缺值处理NaN按前后点线性插值补上 D.load fillmissing(D.load, linear); D.pv fillmissing(D.pv, linear); D.wind fillmissing(D.wind, linear); % 设备参数集中放在para结构体里调参只改这一处 D.para struct(... PG_min, 0.5, ... % 柴油机最小技术出力 MW PG_max, 8.0, ... % 柴油机额定出力 MW ramp, 2.0, ... % 爬坡速率 MW/h ESS_cap, 4.0, ... % 储能额定容量 MWh ESS_Pmax, 1.5, ... % 储能最大充/放电功率 MW eta_c, 0.95, ... % 充电效率 eta_d, 0.92, ... % 放电效率 SOC_min, 0.2, ... % 荷电状态下限 SOC_max, 0.9, ... % 荷电状态上限 a, 0.02, ... % 燃料成本二次项系数 b, 0.5, ... % 燃料成本一次项系数 c, 5.0); % 燃料成本常数项 % 初始荷电状态固定为50%目标函数里循环递推时会用到 D.SOC0 0.5; end几个细节值得注意。fillmissing的linear模式要求被插值向量里至少有两个非NaN相邻点如果数据缺位太多它照样会报错此时要回去检查原始CSV缺了多少行而不是在代码层硬撑。参数集中放在para里也是刻意为之——竞赛后期经常要试好几组储能容量和柴油机上限散落在各函数里的常量每次都要满目录搜索集中定义省掉大量找参数的难堪时间。单位问题也建议在read_data这块统一掉。比如data里给的是kW而目标函数希望用MW计算那读完直接除以1000别在下游每个函数里乘来乘去否则最后调度结果量级一定出错。单位错误是最隐蔽的因为曲线形状看起来完全合理只有成本数字离谱到没办法解释。3.2 目标函数与约束检查把调度模型写进一段代码objective_fun是整套程序的数学核心。它接收一个48维候选解x先拆成PG和PB再算功率平衡、燃料消耗、电网购售成本和约束违反量最后返回一个带罚项的标量fval。PSO就是靠不断尝试不同x、比较fval大小完成寻优。function fval objective_fun(x, D) % x : 48维列向量前24维为柴油机出力PG后24维为储能功率PB % fval: 带罚函数的总成本PSO用它评价每个粒子的优劣 PG x(1:24); PB x(25:48); % 净负荷需要从电网吸收的功率正为购电、负为反送 Pnet D.load - D.pv - D.wind - PB; % 购电成本 售电收入售电按8折结算 buy_cost sum(D.price .* max(Pnet, 0)); sell_inc sum(0.8 * D.price .* max(-Pnet, 0)); grid_cost buy_cost - sell_inc; % 柴油机燃料成本二次函数 a*P^2 b*P c fuel_cost sum(D.para.a .* PG.^2 D.para.b .* PG D.para.c); % 储能SOC递推充电与放电使用不同效率 SOC zeros(25, 1); SOC(1) D.SOC0; for t 1:24 if PB(t) 0 SOC(t1) SOC(t) D.para.eta_c * PB(t) / D.para.ESS_cap; else SOC(t1) SOC(t) - (-PB(t)) / (D.para.eta_d * D.para.ESS_cap); end end % 约束违反量按“超出多少记多少”的方式累积 viol 0; viol viol sum(max(PG - D.para.PG_max, 0)); % 机组超上限 viol viol sum(max(D.para.PG_min - PG, 0)); % 低于下限 viol viol sum(max(abs(diff(PG)) - D.para.ramp, 0)); % 爬坡越限 viol viol sum(max(abs(PB) - D.para.ESS_Pmax, 0)); % 储能功率越限 viol viol sum(max(SOC(1:24) - D.para.SOC_max, 0)); % SOC超上界 viol viol sum(max(D.para.SOC_min - SOC(1:24), 0)); % SOC低于下界 % 罚函数系数让约束违反量在总目标里占据不可忽视的比重 M 1000; fval grid_cost fuel_cost M * viol; end罚函数是这个函数里最值得琢磨的设定。目标成本一般是几百到几千的量级M取1000之后任何一小时的SOC越界都会被放大成几百的目标增量粒子只要碰到不可行解fval就会明显跳高算法就会自然躲开这些区域。这里有个和平时代码习惯不太一样的地方SOC递推用for循环而不是向量化。因为SOC(t1)依赖SOC(t)严格说是时序递推MATLAB数组运算处理不了这种依赖硬向量化反而容易错。一段清楚明白的for循环在竞赛场景远比“看起来高级”的矢量化稳妥。如果想让SOC在一天结束时回到初始值可以在viol上再加abs(SOC(25) - D.SOC0)效果是引导算法最终给出“储能日内复用”的结果。这个小改动在题目要求“储能可持续运行”时非常管用。3.3 PSO求解器收敛与否看这几个参数pso_solve实现的是带线性递减惯性权重的标准粒子群。核心是每个粒子按“自己的历史最优位置”和“全群历史最优位置”两个方向飞迭代次数越大惯性权重越小搜索从全局探索逐渐过渡到局部精细搜索。function [x_best, f_best, curve] pso_solve(obj_fun, dim, lb, ub, opt) % obj_fun: 目标函数句柄 % dim : 决策变量维数 % lb, ub : 变量下界 / 上界向量 % opt : 种群规模、最大迭代、学习因子、惯性权重起止值 swarm opt.swarm_size; maxit opt.max_iter; c1 opt.c1; c2 opt.c2; w0 opt.w0; w1 opt.w1; % 初始化粒子位置和速度均匀分布 小扰动速度 x rand(swarm, dim) .* (ub - lb) lb; v 0.1 * randn(swarm, dim) .* (ub - lb); pbest x; pbest_f arrayfun((i) obj_fun(x(i,:)), 1:swarm); [gbest_f, gidx] min(pbest_f); gbest x(gidx, :); curve zeros(maxit, 1); for it 1:maxit w w0 - (w0 - w1) * it / maxit; % 惯性权重线性递减 for i 1:swarm r1 rand(1, dim); r2 rand(1, dim); v(i,:) w * v(i,:) c1*r1.*(pbest(i,:) - x(i,:)) c2*r2.*(gbest - x(i,:)); x(i,:) x(i,:) v(i,:); % 越界处理直接钳回边界 x(i,:) min(max(x(i,:), lb), ub); fi obj_fun(x(i,:)); if fi pbest_f(i) pbest_f(i) fi; pbest(i,:) x(i,:); if fi gbest_f gbest_f fi; gbest x(i,:); end end end curve(it) gbest_f; end x_best gbest; f_best gbest_f; end惯性权重w00.9、w10.4是一套在竞赛题里很常见的取值范围前几十代粒子飞行速度快有能力跳出局部区域后期w降到0.4左右飞行步长变小在最优附近慢慢精修。c1和c2都取2.0是经典配置含义是“个体经验”和“群体经验”同等重要。越界处理用了“直接钳回边界”而不是把出界粒子随机重置进搜索空间。这条经验来自跑题过程中踩过的坑如果出界后随机重定位粒子瞬间跳出原来的搜索区域爬坡约束这种前后关联的约束被打散收敛曲线会一直抖。钳位之后粒子还被压在那片区域里约束之间的关联性不被破坏收敛稳定得多。真要优化性能可以把内层for循环改成矩阵批量计算但48维、100个粒子、500次迭代个人电脑上也就几分钟。比赛场景里时间多半花在改模型而不是改算法效率上所以这份代码保留逐粒子写法方便对照公式检查。3.4 结果可视化让调度曲线直接能放进论文plot_schedule把最优解从向量里拆出来画成常见两联图上面是负荷、新能源出力、柴油出力、储能放电功率的24小时曲线下面是SOC变化曲线。function plot_schedule(x_best, D, fname) % 绘制日前调度结果并保存为png文件名由调用方传入 PG x_best(1:24); PB x_best(25:48); t 1:24; % 先递推一遍SOC绘图用 SOC zeros(25, 1); SOC(1) D.SOC0; for h 1:24 if PB(h) 0 SOC(h1) SOC(h) D.para.eta_c * PB(h) / D.para.ESS_cap; else SOC(h1) SOC(h) - (-PB(h)) / (D.para.eta_d * D.para.ESS_cap); end end figure(Position, [100 100 900 480]); subplot(2,1,1); plot(t, D.load, k-o, LineWidth, 1.2); hold on; area(t, D.pv D.wind, FaceAlpha, 0.3, EdgeColor, g); stairs(t, PG, r-s, LineWidth, 1.2); stairs(t, -PB, b-d, LineWidth, 1.2); % 放电为负向上显示 legend(负荷, 光伏风电, 柴油出力, 储能放/充, Location, best); xlabel(小时); ylabel(功率/MW); grid on; subplot(2,1,2); stairs(t, SOC(1:24), m-^); yline(D.para.SOC_max, r--); yline(D.para.SOC_min, r--); ylabel(SOC); xlabel(小时); grid on; ylim([0.1 1.0]); saveas(gcf, [../result/, fname]); end出图有几个竞赛常用讲究legend放在不挡曲线的位置网格线开着单位写清楚。论文里最常被问到的“柴油机爬坡”和“储能SOC”恰好是这两张图能直接回答的问题所以图的作用不只是好看还是评阅时的证据。保存路径固定写到../result下。若把这段换成Python版输出要注意matplotlib默认分辨率可能偏小导出前手动设置尺寸和dpi能省去很多排版返工。4. 运行与避坑从解压到收敛这条路值得留一份心4.1 运行顺序与最小验证拿到压缩包后别急着改参数先按原样跑一轮。打开MATLAB把当前路径切到code目录直接运行main.m。机器配置正常时100个粒子、500次迭代大概几十秒到两三分钟。跑完后result目录下会出现schedule_result.png再回命令行看f_best是否和参考值同量级。所谓“最小验证”就是在改任何模型之前先确认原始程序给出的调度结果满足三条直觉负荷大的时段柴油出力不应停在0储能大多在电价低谷充电、尖峰放电SOC始终位于上下限内。这三条任意一条不满足都不要继续往下改先回代码里查约束或数据读取。如果没有参考值我一般会把f_best和手工计算做粗略对照。比如取第10小时把该小时负荷、新能源、柴油出力、储能功率代进目标函数手工算出这一小时成本再和程序内部计算值比一下。这种核对方法很土但对付“单位弄错”“符号反了”基本一抓一个准。4.2 参数调整优先动哪些值PSO对参数不敏感不代表参数随便设就能跑出好结果。实际比赛后期我一般只动四个地方参数默认值调整方向看什么指标M罚函数系数1000增大到5000约束违反是否清零w0惯性权重初值0.9提高到0.95收敛曲线前段是否跳得更开swarm_size种群规模100减少到60单次运行时间max_iter迭代次数500增加到800曲线尾部是否变平经验是先调M确保可行域约束被尊重再调max_iter确保收敛曲线尾部平了最后才动c1、c2和w。上来就调学习因子基本等于对着收敛曲线抓玄学浪费时间。M太小最常见的现象是SOC越界但成本很低粒子觉得“越界也划算”结果里SOC频繁贴到0.1以下。把M调到5000甚至10000约束违反量稍微放大就足以改变粒子决策。M太大的反效果是目标函数变得很陡粒子在可行域边缘震荡收敛曲线像锯齿这时把M降回两三千。4.3 高频问题排查五个问题对应的现象、原因与解决问题一运行read_data时报错“使用readmatrix时出错”。原因基本是MATLAB版本低于R2019a老版本没有readmatrix函数或者CSV第一行带了表头让数据解析失败。解决是先把CSV表头删掉再改用csvread或importdata读取两个办法都能绕开版本兼容问题。问题二程序能跑但负荷曲线形状完全走样比如午间负荷比晚间还高。常见原因是数据单位是kW而代码按MW算或者CSV排列顺序不是0点到23点而是从某个峰谷时段开始。解决方法是先用文本编辑器瞄一眼CSV确认数值量级再交给read_data再把D.load打印到命令行对照24小时曲线基本一眼能看出问题。这算是我最想吃的一颗后悔药因为检查成本太低翻车代价却太大。问题三SOC曲线长时间贴在上边界或下边界。这通常是储能容量ESS_cap与功率ESS_Pmax的比值不协调或者SOC上下限给得太窄。解决是把SOC_min/SOC_max放宽到0.1~0.95同时检查充放电效率是否颠倒——效率写反会让SOC在每次循环里越存越多或越放越少出现违反能量守恒的奇异曲线。问题四收敛曲线前段下降很快后段完全不动且结果明显不是全局最优。这是标准的早熟收敛。解决是在pso_solve里加入粒子随机重初始化逻辑比如每迭代100代挑出10%粒子重新随机撒在搜索空间里让它们去探索新区域通常能跳出局部极值。问题五拿fmincon换掉PSO后直接报NaN。原因是没有光滑化处理的目标函数里有max(Pnet,0)这类不可导项SQP算法的一阶导数和Hessian都算不对。解决一是把max改成softplus近似二是干脆继续用PSO这类不依赖梯度的算法。这也是案例程序坚持用粒子群的底层原因——不追求数学上最漂亮但保证你在赛场上拿得到结果。提示以上五条是跑这类调度案例程序时最常遇到的坑。前三个问题几乎都在数据层后两个才是算法层。比赛时间紧先做数据核对再碰算法参数能省出一整个晚上。5. 二次开发把案例程序改成你自己的调度模型5.1 换数据与加决策变量换数据是把程序私有化最基础的一步。把load24.csv、pv24.csv等四个CSV替换成自己从题目里提取的24小时序列保持列数和单位不变即可。如果题目给定的数据粒度是15分钟那就先把96点聚合成24点或者把决策变量维度从24改成96同时把main.m里的dim、lb、ub全部跟着改目标函数里的时间循环也跟着变。这里最容易遗漏的是爬坡约束15分钟粒度的爬坡限值一般是小时级的四分之一换算不对会得到特别激进的出力曲线。加决策变量是另一个高频需求。比如题目要求考虑可中断负荷可以在现有48维后面追加24维负荷削减比例取值范围0到1% main.m中把dim从48改成72 lb [ones(24,1)*D.para.PG_min; ones(24,1)*(-D.para.ESS_Pmax); zeros(24,1)]; ub [ones(24,1)*D.para.PG_max; ones(24,1)*D.para.ESS_Pmax; ones(24,1)]; % objective_fun中把净负荷的公式改成 Pcut x(49:72); Pnet D.load .* (1 - Pcut) - D.pv - D.wind - PB; % 再在目标函数里补一项失负荷成本 cut_cost sum(D.para.cut_price .* Pcut .* D.load);这种增量式改法风险最小因为原有48维的物理含义完全不动只在新维度上做文章。每加一组变量都建议跑一次原始数据对比加变量前后的成本差是否合理防止新约束把旧结果带偏。5.2 从PSO到整数规划当题目出现0/1变量如果题目进一步要求柴油机启停状态也就是引入0/1变量PSO就不太顺手了。粒子群处理0/1变量通常要先做二进制映射或阈值判断收敛速度明显变慢评阅时解释起来也费劲。此时更稳的办法是换YALMIP加求解器把问题建模成混合整数线性规划% 示意用YALMIP声明24个二元变量代表机组启停 u binvar(24, 1); PG sdpvar(24, 1); Constraints [PG_min .* u PG PG_max .* u, ... ]; % 启停与出力绑定 Objective sum(D.para.a .* PG.^2 D.para.b .* PG D.para.c) ...; optimize(Constraints, Objective);注意目标函数里还有二次项严格说会变成MIQP竞赛中常见做法是先把燃料成本分段线性化再转成MILP求解。这个转换思路建议保留在论文里它是评阅时体现建模功力的关键点。从那以后我每次拿到别人写的调度程序都强制先按原始参数把原数据跑通一遍把目标函数里每一项手工验一遍再动模型换数据只换data目录改模型只动objective_fun和main里的维度声明连result目录都不乱放东西。这个习惯让二次开发从“拆黑匣子”变成“改自己代码”踩坑概率低一大截。希望这种工作方法也能帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/23 17:49:57 · 阅读 1.2k 电工杯B题微电网经济调度MATLAB粒子群求解案例
简介面向2021年电工杯数学建模B题这份资料将赛题思路解析与案例程序打包成套适合备赛选手、指导教师以及所有希望提升数学建模与电气工程应用能力的学习者。资源共包含70个文件压缩后约22.71MB其中37个mat数据文件提供仿真实测数据12个m脚本承载模型求解与绘图代码另有tif/png可视化图片、xlsx中间表格、pdf/doc完整说明文档等便于按数据、算法、结果、报告四类快速定位。目前已有3907人学习浏览。内容上既有对B题背景的深入分析和模型构建思路也有可直接运行调试的Matlab/Python案例程序初学者可借助注释代码理解从实际问题到数学模型的转化进阶者则可参考其算法设计并在电力系统优化、故障诊断等场景中复用。通过研读这些程序与数据读者能够较全面地掌握一类工程建模题的完整解法为后续参赛或科研积累可迁移的实操经验。1. 一份电工杯B题案例程序专治赛期“优化不收敛”2021年电工杯数学建模B题里的微电网经济调度是那种拿到题能列式子、跑起来却反复翻车的典型优化问题。负荷、光伏、风电、分时电价四类数据一进来变量维数、单位换算、爬坡约束和储能荷电状态会同时较劲很多队伍卡在“花了半天把约束写全粒子群一跑却不收敛”。这套案例程序把整个任务拆成数据读取、目标函数、粒子群迭代、出图四个独立模块用MATLAB在几分钟内给出一组可行的24小时调度方案。下面从程序结构开始拆讲清楚怎么跑、怎么改。2. 程序框架与解题思路从赛题描述到48维决策向量2.1 赛题要做的事微电网日前经济调度题目背景就是一个装了柴油发电机、储能、光伏和风电的并网型微电网24小时负荷曲线已知光伏和风电出力按典型日给出分时购电价也定好了。要做的决定是每个小时柴油机发多少、储能充还是放、从电网买多少或卖多少最后让一整天的总运行成本最低。把问题写成数学形式是这样决策变量PG(t)、PB(t)t1,…,24 目标min Σ_t [ a·PG(t)² b·PG(t) c π_buy(t)·max(Pnet(t),0) - 0.8·π_buy(t)·max(-Pnet(t),0) ] 其中 Pnet(t) Pload(t) - Ppv(t) - Pwind(t) - PB(t) 约束PG_min ≤ PG(t) ≤ PG_max |PG(t)-PG(t-1)| ≤ ramp |PB(t)| ≤ ESS_Pmax SOC_min ≤ SOC(t) ≤ SOC_max决策向量拼接成48维前24维是柴油机每小时出力后24维是储能功率。这里有个贯穿全套代码的符号约定——PB大于0表示储能充电等于从负荷侧吸收功率PB小于0表示放电向负荷供电。买电和卖电的价格不等卖电按购电价的80%计算贴近分布式光伏“自发自用、余电上网”的场景。这个模型里有几个容易写错的地方。第一SOC递推不是线性公式充电效率和放电效率不一样放电时要做除法而不是乘法第二爬坡约束是相邻小时之间的差分第1小时没有历史值可以不约束或单独处理第三如果题目要求储能日内复用即每天结束SOC回到初始值那就得在约束或目标里补上“两端守恒”的惩罚项。从竞赛角度看这道题用线性规划也能解Pnet、PG、PB都是连续变量目标函数里带max的部分可以拆成购、售两个非负变量线性化。案例程序之所以选粒子群一方面考虑到不少参赛队对线性规划建模不熟另一方面PSO写起来直观、调参空间大评阅时也能清晰看到模型思考过程不至于被一套黑盒求解器盖住工作量。2.2 解压后的目录结构数据、代码、文档各管一摊整个压缩包解开之后是下面这个布局B_case/ ├── data/ │ ├── load24.csv │ ├── pv24.csv │ ├── wind24.csv │ ├── price24.csv │ └── ESS_param.xlsx ├── code/ │ ├── main.m │ ├── read_data.m │ ├── objective_fun.m │ ├── pso_solve.m │ └── plot_schedule.m ├── result/ │ └── 运行后生成的图和调度数据 └── doc/ └── B题题目说明.pdfdata目录放原始数据code目录放全部脚本result目录在程序运行后自动生成输出图。这类结构在竞赛资源包里很常见核心思想是让数据、逻辑和输出分离——比赛后期换一组数据跑批只改data目录里的CSV代码一行不用动。程序执行顺序是main.m → read_data.m → objective_fun.m → pso_solve.m → plot_schedule.m。read_data返回带全部参数的结构体Dobjective_fun被pso_solve作为函数句柄反复调用最后统一出图。注意所有code里的脚本都用../data/这种相对路径引用原始数据所以别把某个子目录单独拷到别处运行哪怕只是换了一级目录名都可能直接报“找不到文件”。从这个意义上讲解压后保持目录结构原样就是第一个运行前提。main.m的主流程值得先看整个调用关系都在里面% main.m: 一键运行微电网日前经济调度 clear; clc; close all; D read_data(); % 读取24小时基础数据 dim 24 * 2; % 决策变量维数24h柴油出力 24h储能功率 lb [ones(24,1)*D.para.PG_min; ones(24,1)*(-D.para.ESS_Pmax)]; % 下限储能最小功率为负 ub [ones(24,1)*D.para.PG_max; ones(24,1)*D.para.ESS_Pmax]; % 上限储能最大功率为正 opt struct(swarm_size, 100, ... % 种群规模 max_iter, 500, ... % 最大迭代次数 c1, 2.0, c2, 2.0, ... % 个体学习因子 / 全局学习因子 w0, 0.9, w1, 0.4); % 惯性权重的起止值 [x_best, f_best, curve] pso_solve((x) objective_fun(x, D), ... dim, lb, ub, opt); plot_schedule(x_best, D, schedule_result.png); disp([最优成本 , num2str(f_best)]);这个脚本把变量边界直接写在下界/上界向量里而不是塞进目标函数做硬编码这是调参时最省事的写法。后面想加第三类决策变量比如可中断负荷比例只需要把dim改成24*3再在lb和ub里补一段0到1的边界目标函数同步读取就行。2.3 技术栈取舍为什么这里用MATLAB而不是Python很多参赛队第一反应是用Python写赛后复盘也觉得Python库多、更通用但这个包的实现选的是MATLAB原因很实际。一是MATLAB的矩阵语法和这个题目完美匹配。PG、PB、Pload这些全是24维向量mean、max、diff这类运算直接作用在数组上目标函数十几行就能写完。二是绘图方便plot、stairs、area几行就能出一张符合论文格式的图比赛修改比Python里调matplotlib省心。三是评阅老师对MATLAB这套语法太熟悉读代码成本低。四是比赛期间网络不一定稳定Python环境哪怕只是装个pandas和matplotlib也可能碰上镜像源问题MATLAB解压即用省掉一个潜在变量。如果包里有一份Python脚本我一般建议放在“数据初检”这个角色上比如提交前用pandas快速扫一遍CSV里有没有空值和数量级异常再回MATLAB正式跑。正式求解和出图别跨语言切换不然调试时来回换环境非常痛苦。3. 核心代码拆解从CSV到调度曲线的一条链路3.1 数据读取与预处理read_data.mread_data负责把五个文件读进来顺便做一层轻量清洗。CSV原始数据都是24行一列纯数值无表头readmatrix直接读就能拿到24x1向量。若文件带表头或单位注释行readmatrix在R2019a之后会自动把表头识别成字符串行跳过这点在换自己数据时要注意。function D read_data() % read_data 读取微电网24小时数据返回一个结构体D % D.load / D.pv / D.wind / D.price 均为24x1列向量 % D.para 保存柴油机与储能参数 base ../data/; D.load readmatrix([base, load24.csv]); % 24x1负荷功率 MW D.pv readmatrix([base, pv24.csv]); % 24x1光伏出力 MW D.wind readmatrix([base, wind24.csv]); % 24x1风电出力 MW D.price readmatrix([base, price24.csv]); % 24x1购电价 元/MWh % 缺值处理NaN按前后点线性插值补上 D.load fillmissing(D.load, linear); D.pv fillmissing(D.pv, linear); D.wind fillmissing(D.wind, linear); % 设备参数集中放在para结构体里调参只改这一处 D.para struct(... PG_min, 0.5, ... % 柴油机最小技术出力 MW PG_max, 8.0, ... % 柴油机额定出力 MW ramp, 2.0, ... % 爬坡速率 MW/h ESS_cap, 4.0, ... % 储能额定容量 MWh ESS_Pmax, 1.5, ... % 储能最大充/放电功率 MW eta_c, 0.95, ... % 充电效率 eta_d, 0.92, ... % 放电效率 SOC_min, 0.2, ... % 荷电状态下限 SOC_max, 0.9, ... % 荷电状态上限 a, 0.02, ... % 燃料成本二次项系数 b, 0.5, ... % 燃料成本一次项系数 c, 5.0); % 燃料成本常数项 % 初始荷电状态固定为50%目标函数里循环递推时会用到 D.SOC0 0.5; end几个细节值得注意。fillmissing的linear模式要求被插值向量里至少有两个非NaN相邻点如果数据缺位太多它照样会报错此时要回去检查原始CSV缺了多少行而不是在代码层硬撑。参数集中放在para里也是刻意为之——竞赛后期经常要试好几组储能容量和柴油机上限散落在各函数里的常量每次都要满目录搜索集中定义省掉大量找参数的难堪时间。单位问题也建议在read_data这块统一掉。比如data里给的是kW而目标函数希望用MW计算那读完直接除以1000别在下游每个函数里乘来乘去否则最后调度结果量级一定出错。单位错误是最隐蔽的因为曲线形状看起来完全合理只有成本数字离谱到没办法解释。3.2 目标函数与约束检查把调度模型写进一段代码objective_fun是整套程序的数学核心。它接收一个48维候选解x先拆成PG和PB再算功率平衡、燃料消耗、电网购售成本和约束违反量最后返回一个带罚项的标量fval。PSO就是靠不断尝试不同x、比较fval大小完成寻优。function fval objective_fun(x, D) % x : 48维列向量前24维为柴油机出力PG后24维为储能功率PB % fval: 带罚函数的总成本PSO用它评价每个粒子的优劣 PG x(1:24); PB x(25:48); % 净负荷需要从电网吸收的功率正为购电、负为反送 Pnet D.load - D.pv - D.wind - PB; % 购电成本 售电收入售电按8折结算 buy_cost sum(D.price .* max(Pnet, 0)); sell_inc sum(0.8 * D.price .* max(-Pnet, 0)); grid_cost buy_cost - sell_inc; % 柴油机燃料成本二次函数 a*P^2 b*P c fuel_cost sum(D.para.a .* PG.^2 D.para.b .* PG D.para.c); % 储能SOC递推充电与放电使用不同效率 SOC zeros(25, 1); SOC(1) D.SOC0; for t 1:24 if PB(t) 0 SOC(t1) SOC(t) D.para.eta_c * PB(t) / D.para.ESS_cap; else SOC(t1) SOC(t) - (-PB(t)) / (D.para.eta_d * D.para.ESS_cap); end end % 约束违反量按“超出多少记多少”的方式累积 viol 0; viol viol sum(max(PG - D.para.PG_max, 0)); % 机组超上限 viol viol sum(max(D.para.PG_min - PG, 0)); % 低于下限 viol viol sum(max(abs(diff(PG)) - D.para.ramp, 0)); % 爬坡越限 viol viol sum(max(abs(PB) - D.para.ESS_Pmax, 0)); % 储能功率越限 viol viol sum(max(SOC(1:24) - D.para.SOC_max, 0)); % SOC超上界 viol viol sum(max(D.para.SOC_min - SOC(1:24), 0)); % SOC低于下界 % 罚函数系数让约束违反量在总目标里占据不可忽视的比重 M 1000; fval grid_cost fuel_cost M * viol; end罚函数是这个函数里最值得琢磨的设定。目标成本一般是几百到几千的量级M取1000之后任何一小时的SOC越界都会被放大成几百的目标增量粒子只要碰到不可行解fval就会明显跳高算法就会自然躲开这些区域。这里有个和平时代码习惯不太一样的地方SOC递推用for循环而不是向量化。因为SOC(t1)依赖SOC(t)严格说是时序递推MATLAB数组运算处理不了这种依赖硬向量化反而容易错。一段清楚明白的for循环在竞赛场景远比“看起来高级”的矢量化稳妥。如果想让SOC在一天结束时回到初始值可以在viol上再加abs(SOC(25) - D.SOC0)效果是引导算法最终给出“储能日内复用”的结果。这个小改动在题目要求“储能可持续运行”时非常管用。3.3 PSO求解器收敛与否看这几个参数pso_solve实现的是带线性递减惯性权重的标准粒子群。核心是每个粒子按“自己的历史最优位置”和“全群历史最优位置”两个方向飞迭代次数越大惯性权重越小搜索从全局探索逐渐过渡到局部精细搜索。function [x_best, f_best, curve] pso_solve(obj_fun, dim, lb, ub, opt) % obj_fun: 目标函数句柄 % dim : 决策变量维数 % lb, ub : 变量下界 / 上界向量 % opt : 种群规模、最大迭代、学习因子、惯性权重起止值 swarm opt.swarm_size; maxit opt.max_iter; c1 opt.c1; c2 opt.c2; w0 opt.w0; w1 opt.w1; % 初始化粒子位置和速度均匀分布 小扰动速度 x rand(swarm, dim) .* (ub - lb) lb; v 0.1 * randn(swarm, dim) .* (ub - lb); pbest x; pbest_f arrayfun((i) obj_fun(x(i,:)), 1:swarm); [gbest_f, gidx] min(pbest_f); gbest x(gidx, :); curve zeros(maxit, 1); for it 1:maxit w w0 - (w0 - w1) * it / maxit; % 惯性权重线性递减 for i 1:swarm r1 rand(1, dim); r2 rand(1, dim); v(i,:) w * v(i,:) c1*r1.*(pbest(i,:) - x(i,:)) c2*r2.*(gbest - x(i,:)); x(i,:) x(i,:) v(i,:); % 越界处理直接钳回边界 x(i,:) min(max(x(i,:), lb), ub); fi obj_fun(x(i,:)); if fi pbest_f(i) pbest_f(i) fi; pbest(i,:) x(i,:); if fi gbest_f gbest_f fi; gbest x(i,:); end end end curve(it) gbest_f; end x_best gbest; f_best gbest_f; end惯性权重w00.9、w10.4是一套在竞赛题里很常见的取值范围前几十代粒子飞行速度快有能力跳出局部区域后期w降到0.4左右飞行步长变小在最优附近慢慢精修。c1和c2都取2.0是经典配置含义是“个体经验”和“群体经验”同等重要。越界处理用了“直接钳回边界”而不是把出界粒子随机重置进搜索空间。这条经验来自跑题过程中踩过的坑如果出界后随机重定位粒子瞬间跳出原来的搜索区域爬坡约束这种前后关联的约束被打散收敛曲线会一直抖。钳位之后粒子还被压在那片区域里约束之间的关联性不被破坏收敛稳定得多。真要优化性能可以把内层for循环改成矩阵批量计算但48维、100个粒子、500次迭代个人电脑上也就几分钟。比赛场景里时间多半花在改模型而不是改算法效率上所以这份代码保留逐粒子写法方便对照公式检查。3.4 结果可视化让调度曲线直接能放进论文plot_schedule把最优解从向量里拆出来画成常见两联图上面是负荷、新能源出力、柴油出力、储能放电功率的24小时曲线下面是SOC变化曲线。function plot_schedule(x_best, D, fname) % 绘制日前调度结果并保存为png文件名由调用方传入 PG x_best(1:24); PB x_best(25:48); t 1:24; % 先递推一遍SOC绘图用 SOC zeros(25, 1); SOC(1) D.SOC0; for h 1:24 if PB(h) 0 SOC(h1) SOC(h) D.para.eta_c * PB(h) / D.para.ESS_cap; else SOC(h1) SOC(h) - (-PB(h)) / (D.para.eta_d * D.para.ESS_cap); end end figure(Position, [100 100 900 480]); subplot(2,1,1); plot(t, D.load, k-o, LineWidth, 1.2); hold on; area(t, D.pv D.wind, FaceAlpha, 0.3, EdgeColor, g); stairs(t, PG, r-s, LineWidth, 1.2); stairs(t, -PB, b-d, LineWidth, 1.2); % 放电为负向上显示 legend(负荷, 光伏风电, 柴油出力, 储能放/充, Location, best); xlabel(小时); ylabel(功率/MW); grid on; subplot(2,1,2); stairs(t, SOC(1:24), m-^); yline(D.para.SOC_max, r--); yline(D.para.SOC_min, r--); ylabel(SOC); xlabel(小时); grid on; ylim([0.1 1.0]); saveas(gcf, [../result/, fname]); end出图有几个竞赛常用讲究legend放在不挡曲线的位置网格线开着单位写清楚。论文里最常被问到的“柴油机爬坡”和“储能SOC”恰好是这两张图能直接回答的问题所以图的作用不只是好看还是评阅时的证据。保存路径固定写到../result下。若把这段换成Python版输出要注意matplotlib默认分辨率可能偏小导出前手动设置尺寸和dpi能省去很多排版返工。4. 运行与避坑从解压到收敛这条路值得留一份心4.1 运行顺序与最小验证拿到压缩包后别急着改参数先按原样跑一轮。打开MATLAB把当前路径切到code目录直接运行main.m。机器配置正常时100个粒子、500次迭代大概几十秒到两三分钟。跑完后result目录下会出现schedule_result.png再回命令行看f_best是否和参考值同量级。所谓“最小验证”就是在改任何模型之前先确认原始程序给出的调度结果满足三条直觉负荷大的时段柴油出力不应停在0储能大多在电价低谷充电、尖峰放电SOC始终位于上下限内。这三条任意一条不满足都不要继续往下改先回代码里查约束或数据读取。如果没有参考值我一般会把f_best和手工计算做粗略对照。比如取第10小时把该小时负荷、新能源、柴油出力、储能功率代进目标函数手工算出这一小时成本再和程序内部计算值比一下。这种核对方法很土但对付“单位弄错”“符号反了”基本一抓一个准。4.2 参数调整优先动哪些值PSO对参数不敏感不代表参数随便设就能跑出好结果。实际比赛后期我一般只动四个地方参数默认值调整方向看什么指标M罚函数系数1000增大到5000约束违反是否清零w0惯性权重初值0.9提高到0.95收敛曲线前段是否跳得更开swarm_size种群规模100减少到60单次运行时间max_iter迭代次数500增加到800曲线尾部是否变平经验是先调M确保可行域约束被尊重再调max_iter确保收敛曲线尾部平了最后才动c1、c2和w。上来就调学习因子基本等于对着收敛曲线抓玄学浪费时间。M太小最常见的现象是SOC越界但成本很低粒子觉得“越界也划算”结果里SOC频繁贴到0.1以下。把M调到5000甚至10000约束违反量稍微放大就足以改变粒子决策。M太大的反效果是目标函数变得很陡粒子在可行域边缘震荡收敛曲线像锯齿这时把M降回两三千。4.3 高频问题排查五个问题对应的现象、原因与解决问题一运行read_data时报错“使用readmatrix时出错”。原因基本是MATLAB版本低于R2019a老版本没有readmatrix函数或者CSV第一行带了表头让数据解析失败。解决是先把CSV表头删掉再改用csvread或importdata读取两个办法都能绕开版本兼容问题。问题二程序能跑但负荷曲线形状完全走样比如午间负荷比晚间还高。常见原因是数据单位是kW而代码按MW算或者CSV排列顺序不是0点到23点而是从某个峰谷时段开始。解决方法是先用文本编辑器瞄一眼CSV确认数值量级再交给read_data再把D.load打印到命令行对照24小时曲线基本一眼能看出问题。这算是我最想吃的一颗后悔药因为检查成本太低翻车代价却太大。问题三SOC曲线长时间贴在上边界或下边界。这通常是储能容量ESS_cap与功率ESS_Pmax的比值不协调或者SOC上下限给得太窄。解决是把SOC_min/SOC_max放宽到0.1~0.95同时检查充放电效率是否颠倒——效率写反会让SOC在每次循环里越存越多或越放越少出现违反能量守恒的奇异曲线。问题四收敛曲线前段下降很快后段完全不动且结果明显不是全局最优。这是标准的早熟收敛。解决是在pso_solve里加入粒子随机重初始化逻辑比如每迭代100代挑出10%粒子重新随机撒在搜索空间里让它们去探索新区域通常能跳出局部极值。问题五拿fmincon换掉PSO后直接报NaN。原因是没有光滑化处理的目标函数里有max(Pnet,0)这类不可导项SQP算法的一阶导数和Hessian都算不对。解决一是把max改成softplus近似二是干脆继续用PSO这类不依赖梯度的算法。这也是案例程序坚持用粒子群的底层原因——不追求数学上最漂亮但保证你在赛场上拿得到结果。提示以上五条是跑这类调度案例程序时最常遇到的坑。前三个问题几乎都在数据层后两个才是算法层。比赛时间紧先做数据核对再碰算法参数能省出一整个晚上。5. 二次开发把案例程序改成你自己的调度模型5.1 换数据与加决策变量换数据是把程序私有化最基础的一步。把load24.csv、pv24.csv等四个CSV替换成自己从题目里提取的24小时序列保持列数和单位不变即可。如果题目给定的数据粒度是15分钟那就先把96点聚合成24点或者把决策变量维度从24改成96同时把main.m里的dim、lb、ub全部跟着改目标函数里的时间循环也跟着变。这里最容易遗漏的是爬坡约束15分钟粒度的爬坡限值一般是小时级的四分之一换算不对会得到特别激进的出力曲线。加决策变量是另一个高频需求。比如题目要求考虑可中断负荷可以在现有48维后面追加24维负荷削减比例取值范围0到1% main.m中把dim从48改成72 lb [ones(24,1)*D.para.PG_min; ones(24,1)*(-D.para.ESS_Pmax); zeros(24,1)]; ub [ones(24,1)*D.para.PG_max; ones(24,1)*D.para.ESS_Pmax; ones(24,1)]; % objective_fun中把净负荷的公式改成 Pcut x(49:72); Pnet D.load .* (1 - Pcut) - D.pv - D.wind - PB; % 再在目标函数里补一项失负荷成本 cut_cost sum(D.para.cut_price .* Pcut .* D.load);这种增量式改法风险最小因为原有48维的物理含义完全不动只在新维度上做文章。每加一组变量都建议跑一次原始数据对比加变量前后的成本差是否合理防止新约束把旧结果带偏。5.2 从PSO到整数规划当题目出现0/1变量如果题目进一步要求柴油机启停状态也就是引入0/1变量PSO就不太顺手了。粒子群处理0/1变量通常要先做二进制映射或阈值判断收敛速度明显变慢评阅时解释起来也费劲。此时更稳的办法是换YALMIP加求解器把问题建模成混合整数线性规划% 示意用YALMIP声明24个二元变量代表机组启停 u binvar(24, 1); PG sdpvar(24, 1); Constraints [PG_min .* u PG PG_max .* u, ... ]; % 启停与出力绑定 Objective sum(D.para.a .* PG.^2 D.para.b .* PG D.para.c) ...; optimize(Constraints, Objective);注意目标函数里还有二次项严格说会变成MIQP竞赛中常见做法是先把燃料成本分段线性化再转成MILP求解。这个转换思路建议保留在论文里它是评阅时体现建模功力的关键点。从那以后我每次拿到别人写的调度程序都强制先按原始参数把原数据跑通一遍把目标函数里每一项手工验一遍再动模型换数据只换data目录改模型只动objective_fun和main里的维度声明连result目录都不乱放东西。这个习惯让二次开发从“拆黑匣子”变成“改自己代码”踩坑概率低一大截。希望这种工作方法也能帮到你。本文还有配套的精品资源点击获取