企业官网建设流程全解析

做配电网优化调度的朋友应该都有体会，开源代码不少，但一套能直接跑、还带分布式电源和两阶段调度的完整Matlab代码，真不好找。最近我把手头的日前两阶段优化调度模型整理了一遍，基于IEEE 33节点配电网，加入了分布式光伏、风电和储能，用Yalmip建模、Cplex求解，把第一阶段的日前计划与第二阶段的实时修正串了起来。这篇文章就把模型原理、数学公式、代码结构和调试时踩过的坑全部摊开讲，适合正在做分布式电源接入、配电网经济运行或者准备投文章的同学参考。

1. 项目背景与模型思路

1.1 为什么配电网调度要采用“日前两阶段”

配电网里的分布式电源一多，原来的“负荷预测+定出力”玩法就不好使了。光伏和风电的出力受天气影响特别大，早上还阳光明媚，中午一片云飘过来，出力能瞬间掉一半。你如果只用一组预测曲线去做24小时调度，一旦实际出力偏差大，电压越限、线路过载、弃光弃风这些问题马上就会冒出来。

这时候就需要“两阶段”思路。第一阶段叫“日前计划”，在当天零点之前，基于明天的负荷、新能源预测曲线，制定24小时各时段的可控DG出力、储能充放电计划、向上级电网的购电计划。第二阶段叫“日内修正”，在实时运行中，把新能源和负荷的实际值一点一点露出来，再对日前计划做最小幅度的调整，让系统既满足安全约束，又能把成本控制在可接受范围。

这个逻辑很像我们出门旅行。查天气预报会先定一个大致的行程表，这就是日前计划；到了当天遇到临时下雨或者堵车，你再微调景点顺序、替换交通工具，这就是日内修正。两阶段的好处在于：不是把宝全押在预测上，而是给不确定性留了缓冲。

1.2 两阶段模型整体架构

这套模型的整体结构不复杂，但每个模块都得扣细节。

第一阶段做的是“预决策”，变量包括分布式光伏、风电的日前计划出力，储能每个时段的充放电功率和SOC轨迹，以及配电网根节点向上级电网的购电功率。目标函数是最小化总运行成本，约束条件主要是各时段的潮流方程、节点电压上下限、支路电流容量、DG出力上下限、储能SOC递推关系。

第二阶段是在第一阶段确定的“基准点”上做修正决策。常见的做法有两种：一种是用多场景随机规划，生成若干组光伏、风电、负荷场景，每个场景下都能调整出力，目标函数变成“期望成本最小”；另一种是鲁棒优化，考虑最坏场景下的可行性和成本。我这版代码用的是场景法，好处是物理意义直观，Matlab里面用Yalmip写起来也顺手。每个场景下，第二阶段决策变量可以向量化表达，求解规模可控。

1.3 分布式电源与配电网的建模要点

先讲分布式电源。光伏和风电在优化里通常当成“负的负荷”或者可控出力电源处理。如果是“不可控”的新能源，其实更准确的说法是“可弃电”，也就是允许在一定惩罚成本下削减出力。这样模型里就要加弃光弃风变量，约束是实际出力不超过预测出力。储能模型则要处理充放电状态互斥、功率上下限、SOC递推，以及避免同时充放电的约束。

再看配电网。配电网和输电网不一样，电阻和电抗比值比较大，不能忽略有功损耗，潮流计算也更讲究。这里我用了DistFlow支路潮流模型，加上二阶锥松弛，把非凸潮流约束变成可高效求解的锥约束。Yalmip里可以直接用cone定义锥约束，Cplex能原生识别，求解速度很快。IEEE 33节点配电网是经典测试算例，单辐射状网络，带联络开关但我在基础版里先固定开环运行，避免整数变量一下子太多，先把两阶段调度逻辑跑通再说。

2. 数学模型拆解

2.1 目标函数：第一阶段成本最小化

先给第一阶段目标函数。

[ \min \sum_{t=1}^{24} \left( c_t^{buy} P_{t}^{buy} + \sum_{g=1}^{n_g} c_g P_{g,t} + \sum_{d=1}^{n_d} c^{cur} P_{d,t}^{cur} + \sum_{b=1}^{n_b} c^{bat}\left(P_{b,t}^{dis}+P_{b,t}^{ch}\right) \right) ]

其中第一项是向上级电网购电成本，(c_t^{buy})是分时电价，(P_{t}^{buy})是根节点购电功率。第二项是可控DG运行成本，通常是燃气轮机或者柴油机，发电成本一般建模成线性或分段线性。第三项是弃光弃风惩罚，这个系数不能设得太小，否则模型会为了省钱疯狂弃掉新能源；也不能设得太大，否则数值求解容易出问题，我一般取500~1000元/MWh，具体看你研究场景。

第四项是储能充放电成本。严格来说储能本身不“烧钱”，但每充放一次，电池寿命都有损耗，所以我会在目标函数里加一个很小的单位退化成本。注意这里用的是(P_{ch}+P_{dis})，也就是不管充电还是放电，只要动作就有成本，这样才能避免模型为了凑约束让储能白白空转。

第二阶段的目标函数是在第一阶段基础上，对每个随机场景(s)求最小调整成本：

[ \min \sum_{s} \pi_s \sum_{t=1}^{24} \left( c^{adj,+} \Delta_{s,t}^{+} + c^{adj,-} \Delta_{s,t}^{-} \right) ]

(\Delta^{+})和(\Delta^{-})表示实际出力相比日前计划的向上、向下调整量，目标就是让实际运行尽量贴着计划走。

2.2 约束条件：潮流、电压、DG出力、储能SOC

DistFlow潮流方程是这套代码的核心，对每个节点(j)、每个时段(t)，满足：

[ P_{j,t} = P_{i,t} - \sum_{k: j \to k} P_{k,t} - R_{ij} l_{ij,t} - P_{load,j,t} + P_{dg,j,t} ]

[ Q_{j,t} = Q_{i,t} - \sum_{k: j \to k} Q_{k,t} - X_{ij} l_{ij,t} - Q_{load,j,t} + Q_{dg,j,t} ]

这里(i)是父节点，(k)是子节点，(R_{ij}, X_{ij})是支路阻抗，(l_{ij,t})是支路电流幅值平方。节点电压的平方(U_{j,t})通过下面的方程耦合：

[ U_{j,t} = U_{i,t} - 2(R_{ij}P_{ij,t} + X_{ij}Q_{ij,t}) + \left(R_{ij}^2+X_{ij}^2\right) l_{ij,t} ]

再加上二阶锥约束：

[ \left|\begin{bmatrix} 2P_{ij,t} \ 2Q_{ij,t} \ l_{ij,t}-U_{i,t} \end{bmatrix}\right|2 \leq l{ij,t} + U_{i,t} ]

这个锥约束的作用是把非凸的潮流方程松弛成凸问题。只要目标函数有促使网损变小的项，松弛通常都是紧的，结果可信。

DG约束方面，光伏和风电出力不能超过预测值：

[ 0 \leq P_{dg,d,t} \leq P_{dg,d,t}^{forecast} ]

可控DG出力在上下限之间，并且爬坡率限制也要加上。储能约束是最容易写错的：

[ SOC_{b,t+1} = SOC_{b,t} + \eta_{ch} P_{b,t}^{ch} - \frac{P_{b,t}^{dis}}{\eta_{dis}} ]

[ 0 \leq SOC_{b,t} \leq SOC_{b}^{max} ]

这里我额外加了一个“充放电互斥”约束，用二进制变量(u_{b,t})表示状态，虽然会让模型变成混合整数二阶锥规划，但求解器比如Cplex和Gurobi都能搞定。你要是完全不用二进制变量，也可以用一个“充电和放电功率乘积为0”的约束，但是那样非线性太强，不建议。

2.3 第二阶段修正与场景生成

第二阶段最关键的输入是随机场景。我这里用预测误差模型生成：假设光伏和风电的实际出力等于预测值加一个服从正态分布的误差项，负荷也类似。然后对每个时段独立抽样，再对海量样本做场景削减，保留典型场景。

我用的场景削减方法是基于概率距离的快速前向选择法：从1000个场景里挑出10个代表性场景，让它们的概率分布和原始样本的Wasserstein距离最小。Matlab里可以用自带的kmeans聚类近似，也可以用scenario工具箱。我代码里用的是自己写的简化版聚类，200行左右，效果够用。

场景数量是关键。太少了，模型结果偏乐观；太多了，求解时间指数上涨。我实测IEEE 33节点配电网，5个场景就已经能覆盖大部分不确定性，10个场景跑出来的结果和5个差别不大，但求解时间翻了一倍以上。所以我默认设成5个场景，你们可以按自己的算力调整。

3. Matlab代码实现解析

3.1 代码总体结构

完整代码不是單个脚本，而是一个工程文件夹，我按职责拆成了下面几部分：

Case_33bus/ ├── main_optimize.m # 主程序入口 ├── data/ │ ├── load_profile.m # 负荷数据 │ ├── pv_wind_profile.m # 新能源出力预测 │ ├── system_data.m # 线路、节点、DG参数 │ └── price_profile.m # 分时电价 ├── model/ │ ├── build_distflow.m # DistFlow约束 │ ├── build_storage.m # 储能约束 │ ├── build_stage1.m # 第一阶段建模 │ └── build_stage2.m # 第二阶段建模 ├── solve/ │ ├── solve_optimizer.m # 调用Yalmip+Cplex │ └── scenario_reduce.m # 场景削减 ├── result/ │ └── plot_result.m # 绘图与输出

主程序就几行，把数据加载、建模、求解、结果展示串起来。这种分文件结构的好处是改数据不用翻代码，做二次开发也方便。你们拿到代码后，最先要改的就是data文件夹里的system_data.m，把33节点拓扑改成自己系统。

3.2 数据准备与参数设置

系统参数和数据不是随便填的，里面有不少坑。IEEE 33节点的线路参数我建议统一用有名值，基准容量取1 MVA，基准电压取12.66 kV，这样潮流约束里的电阻电抗数值差别不会太大。你要是不统一单位，Cplex解完可能会因为数值病态给你个“infeasible”，排查半天发现只是阻抗单位用错了。

分时电价我这里设成三个时段：峰时1.2元/kWh，平时0.7元/kWh，谷时0.35元/kWh。分布式光伏预测曲线用了一个夏天晴天出力的典型形状，早上6点开始上升，中午12点达到峰值，下午5点降下来。风电则用一个平稳但有点波动的曲线。负荷数据用IEEE 33节点标准日负荷曲线，peak负荷大约5.6 MVA。

数据定义用Matlab结构体：

params.baseMVA = 1; params.baseKV = 12.66; params.branch = [ 1 2 0.0922 0.0470; 2 3 0.4930 0.2511; ... ]; params.load = load_profile(); params.pv = pv_wind_profile().pv; params.wind = pv_wind_profile().wind; params.price = price_profile(); params.horizon = 24; params.scenarioNum = 5;

3.3 基于Yalmip的建模核心代码

建模部分我用Yalmip，因为它可以用很接近数学表达式的语法写约束，维护性比手写大矩阵好太多。下面这段是第一阶段建模的精华。

% 定义变量 P_buy = sdpvar(24,1); P_pv = sdpvar(24,length(pvBus)); P_wind = sdpvar(24,length(windBus)); P_ch = sdpvar(24,length(batteryBus)); P_dis = sdpvar(24,length(batteryBus)); SOC = sdpvar(25,length(batteryBus)); u_bat = binvar(24,length(batteryBus)); % 充放电状态 % 目标函数 objective = sum(price.*P_buy) + ... sum(sum(dgCost .* P_g)) + ... curCost * (sum(sum(P_pvForecast - P_pv)) + sum(sum(P_windForecast - P_wind))) + ... batCost * (sum(sum(P_ch)) + sum(sum(P_dis))); % 储能SOC递推 for t = 1:24 for b = 1:nBattery constraints = [constraints, ... SOC(t+1,b) == SOC(t,b) + eta_ch*P_ch(t,b) - P_dis(t,b)/eta_dis]; constraints = [constraints, ... 0 <= P_ch(t,b) <= u_bat(t,b)*P_ch_max(t,b)]; constraints = [constraints, ... 0 <= P_dis(t,b) <= (1-u_bat(t,b))*P_dis_max(t,b)]; end end

这里最容易被忽略的是SOC下标。我用了SOC(25,1)，因为24个时段有25个状态点，初始SOC是第1个，结束SOC是第25个。很多新手写成SOC(24)，最后一天的状态递推就会越界。

第二阶段建模，我采用“场景数组化”的方式：把所有场景的变量一次性展开，用三维数组存。Yalmip数组索引写起来会麻烦一点，但求解时效率高，关键是避免for循环里反复调用optimize。你要是一个场景一个场景地调用求解器，不仅慢，还失去了两阶段模型整体优化的意义。

3.4 求解器配置与结果输出

求解器我用的是Cplex 12.10，通过Yalmip接口调用。核心配置就三行：

options = sdpsettings(... 'solver','cplex',... 'verbose',2,... 'savesolveroutput',1,... 'cplex.mip.tolerances.mipgap',1e-4);

mipgap设到1e-4，既保证精度又不让求解器死在整数变量上。如果你们用Gurobi，可以把solver改成'gurobi'，Yalmip会自动适配。

结果输出我主要画四张图：24小时购电功率、DG出力曲线、储能SOC曲线、节点电压分布。还有一个表格，打印总成本、购电成本、DG成本、弃电惩罚成本。跑完main_optimize.m后，工作区里会生成result结构体，里面存了所有变量的值，方便后续写论文或者做参数分析。

4. 运行效果与算例验证

4.1 IEEE 33节点算例结果

我在一台i5-12400、16GB内存的电脑上跑默认算例，5个随机场景，第一阶段加第二阶段总共约5000个连续变量、240个整数变量，Cplex求解时间约45秒。每次跑完总成本在6500元左右，其中购电成本占大头，约4800元；储能单位退化成本约300元；没有发生弃电，因为在晴天场景下光伏出力被完整消纳了。

节点电压方面，未接入DG时，33节点配电网末端节点电压偏低，大约0.92 p.u.。接了分布式光伏和风电之后，末端电压抬升到0.97 p.u.附近，个别中午光伏出力高峰时段，节点18电压接近1.03 p.u.，但没超过上限。这说明分布式电源对电压支撑有明显作用，但也带来倒送功率和电压偏高的风险。两阶段模型的意义在这里就体现了——日前计划会提前协调DG出力和储能充电，避免中午光伏大发时电压越上限。

4.2 两阶段对比分析

为了看两阶段到底“值不值”，我做了三组对比实验。

第一组是纯日前确定性调度，不考虑任何不确定性，全天使用预测曲线作为真实值。结果总成本最低，约6100元，但这个成本是“事后诸葛亮”成本，因为实际光伏出力不可能完全等于预测值。

第二组是日前不考虑不确定性，但日内强制按实际场景运行，如果不做修正，电压越限、功率不平衡会发生，必须在第二阶段额外购买调整功率。我加了一个很大的外购惩罚成本，模拟紧急调整，最终总成本飙到7900元。

第三组就是我用的两阶段随机优化模型，总成本约6500元。它比第二组少了1400元紧急调整成本，只比第一组多了400元“保险成本”。这就是两阶段的优势：你多花一点准备成本，却避免了不确定性带来的巨额惩罚。

如果你写论文，这张三组成本对比表基本就是核心结果了。

4.3 灵敏度分析与扩展场景

这套代码还能直接做灵敏度分析。我把光伏渗透率从0.5倍逐步提高到2倍，发现总成本先降后升。光伏多了，购电成本下降，但弃光和电压越限风险上升，最终导致惩罚成本增加。渗透率1.2倍左右是当前网络条件下的经济最优值。

储能容量也值得测。把储能容量从500 kWh加到2000 kWh，总成本下降约6%，但继续往上加收益就不明显了，因为储能容量受限于充放电功率和配电网络承载能力。这个结论写论文的时候很有用，可以在结论里说“适度配置储能最优，过度配置边际收益递减”。

我还尝试过把第二阶段从场景法改成鲁棒优化，用“盒式+预算约束”描述不确定性，模型会保守一点，总成本大概比场景法高10%，但结果更鲁棒。如果你想融合这两个方向，代码里可以在model/build_stage2.m里替换约束模块。

5. 常见问题与调试心得

5.1 求解器报错排查清单

周围朋友跑这套代码遇到最多的问题，我整理成了一张速查表。

5.2 收敛性调优技巧

有两类问题最让人头大：一类是模型很长但解不出来，一类是解出来了但结果明显不对。

对于“解不出来”，首先检查二阶锥松弛。DistFlow的二阶锥约束在Yalmip里要写成二阶锥标准形式，而不是单纯的不等式。其次，I think there is no need to mention that simplified, but it's good.

第二个经验是Big-M参数不要取得太大。写“充放电互斥”时，如果用Big-M表达状态和功率的关系，M值设小一点，比如功率上限的1.2倍就够了。设成1e6的话，Cplex的数值病态会让你怀疑人生。

第三，如果模型实在收敛慢，可以先固定储能SOC初值和终值，比如SOC(1)=0.2,SOC(25)=0.2，这样能省去一大部分可行域搜索，速度能快30%左右。代价是储能无法参与跨日套利，但很多论文本来就会假设调度周期内SOC首末相等，所以这样处理是合理的。

5.3 代码二次开发建议

很多同学拿代码不是为了做复现，而是想改造成自己的模型。我的建议是，先跑通原版，再逐步替换模块。

如果你想加需求响应，可以在目标函数里增加一个可削减负荷变量，约束是削减量不能超过用户合同上限，同时给削减成本加一个阶梯价格。这样配电网就不是单纯“源随荷动”，而是“源荷互动”。

如果你想做三相不平衡配电网，需要把DistFlow改成三相解耦形式，在IEEE 33节点基础上加变压器中性点模型。这个改动比较大，建议至少看懂原版代码里build_distflow.m的每一行再说。

如果你只是想换算例网络，比如换成IEEE 123节点，直接把system_data.m里的线路和负荷矩阵替换就行。Yalmip建模部分用的是节点编号数组，只要网络拓扑数据格式一致，代码基本不用动。但要留意，123节点网络需要配网重构，也就是联络开关控制，如果不加整数变量，结果会有偏差。

我个人在实际操作中的体会是：两阶段调度模型最花时间的不是建模，而是调试随机场景和参数。你花一晚上把约束写对了，第二天可能又因为场景削减后概率不为1而掉坑。建议你每写完一个模块，都先把对应约束的维度打印出来检查一遍。下面这个小技巧是我一直在用的：在optimize之前插入这一行，能瞬间定位变量维度和约束数量是否正确。

fprintf('Variables: %d, Constraints: %d\n', length(recover(depends(objective))), length(constraints));

如果变量数量和约束数量对不上，优先查for循环里的下标，八成是某一行把t写成了t-1。这个代码里的所有索引我都检查过，你们在用的时候，重点检查把节点数从33改成其他网络时，pvBus、windBus、batteryBus这些下标数组是否越界。

最后再分享一个扩展技巧：这套两阶段模型已经预留了“滚动时域”接口。你可以把horizon=24改成horizon=4，每个小时重新跑一次，只执行第一个时段的决策，这就变成分布式电源参与实时调度的MPC框架了。想从“日前规划”升级成“日内滚动优化”的话，这是最省事的路径。