手把手教你用 MATLAB 实现 PSO + CPLEX 双层优化（附完整数学模型与代码架构）

🚀 1. 前言

在综合能源系统（IES）、需求响应（DR）、微电网优化调度以及电力市场交易等研究领域中，我们经常会遇到这样一类“套娃”问题：

• 运营商（上层）：希望通过调整电价、热价或碳价，割一波韭菜……啊不，是最大化自身运行收益。
• 用户端（下层）：面对运营商发布的价格信号，见招拆招，优化自身的设备出力和用能策略，以最小化自身用能成本。

双方存在着明显的主从博弈关系（Leader-Follower）。传统的单层优化对此类问题无能为力，因而双层优化（Bi-level Optimization）应运而生。

在实际工程和论文中，上层决策变量（如动态定价）往往具有非线性、非凸特征，而下层调度（如设备物理约束）通常可以构建为混合整数线性规划（MILP）模型。因此，行业内诞生了一种黄金求解组合：

💡PSO（粒子群算法）+ CPLEX（数学规划求解器）

• 外层（智能算法）：采用 PSO 负责全局搜索最优价格策略。
• 内层（精确求解）：采用 CPLEX 负责秒杀下层的设备最优调度方案。

本文将从主从博弈原理、数学模型构建、MATLAB代码工程设计三个维度，手把手带你调通这个双层优化框架！

💡 2. 双层优化的核心：主从博弈关系

双层优化本质上是一个Stackelberg 博弈模型。上层领导者先出牌，下层跟随者根据上层的牌局做出最优反应，其数学抽象表达如下：

在综合能源系统（IES）的落地场景中，这个框架的职责分工非常明确：

角色	核心主体	决策变量	优化目标
上层 (Leader)	能源运营商	电价、热价、碳价/补贴	（购销差价利润最大化）
下层 (Follower)	终端用户/微网	负荷响应、储能充放、P2G与CCS出力	（综合用能成本最小化）

📐 3. 上层优化模型构建

以 24 小时动态电价制定为例，上层运营商通过控制价格来引导下层负荷。

3.1 决策变量

上层粒子的位置编码即为 24 小时的电价曲线：

3.2 目标函数

运营商的净利润等于销售电费收入减去向大电网购电的成本：

其中，为大电网电价（已知量），和均需要由下层优化求解后返回给上层。

📐 4. 下层优化模型构建

下层用户在接收到上层发布的电价后，开始规划自己的一日用能策略。

4.1 目标函数

用户追求综合运行成本最小化：

各个成本子项的数学逻辑如下：

• 电费成本：
• 碳交易成本：
• 用电舒适度损失（可调负荷惩罚）：

4.2 核心系统约束

为了保证物理系统不崩溃，下层需要严格满足各类等式与不等式约束：

• 电功率平衡约束：

• 储能SOC动态约束：

• P2G（电转气）与 CCS（碳捕集）耦合约束：

🛠️ 5. MATLAB 工程结构设计

在写代码之前，良好的工程目录习惯能让你事半功倍。推荐采用模块化设计，结构如下：

1Project (双层优化项目主目录) 2│ 3├── main.m % 主程序（初始化参数、调用PSO、打印最终画图） 4├── PSO.m % 粒子群算法主循环核心 5├── upper_obj.m % 上层目标函数（负责搭建上下层传递的桥梁） 6├── lower_cplex.m % 下层CPLEX求解器（基于YALMIP工具箱搭建） 7├── plot_result.m % 结果可视化函数 8└── data.mat % 风光负荷基础数据

💻 6. 核心代码手把手实现

6.1 上层适应度计算 (`upper_obj.m`)

这个函数是连接上下层的纽带。它接收 PSO 传过来的当前“价格方案”，丢给下层 CPLEX 算一遍，再计算出运营商的利润。

1function minus_profit = upper_obj(price) 2 % 输入：price - 当前粒子代表的24小时电价策略 3 4 % 调用下层 CPLEX 模型求解 5 result = lower_cplex(price); 6 7 % 提取下层响应后的最优调度结果 8 load_respond = result.load; 9 grid_buy = result.grid; 10 grid_price = result.grid_price; % 大电网已知电价 11 12 % 计算运营商利润 (销售收入 - 购电成本) 13 profit = sum(price .* load_respond) - sum(grid_price .* grid_buy); 14 15 % 由于标准PSO通常寻找最小值，因此将最大化利润取负号 16 minus_profit = -profit; 17end

6.2 下层智能调度模型 (`lower_cplex.m`)

下层基于YALMIP + CPLEX搭建。记得把上层传下来的price当作已知常数处理。

1function result = lower_cplex(price) 2 % 1. 声明决策变量 (24小时出力) 3 Pgrid = sdpvar(24,1); % 购电 4 Pch = sdpvar(24,1); % 充电 5 Pdis = sdpvar(24,1); % 放电 6 SOC = sdpvar(24,1); % 蓄电池容量状态 7 8 % 2. 导入基础数据 (此处简写，实际从data.mat加载) 9 PV = ...; WT = ...; Load = ...; 10 eta_c = 0.95; eta_d = 0.95; 11 12 % 3. 构建目标函数 13 obj = 0; 14 for t = 1:24 15 obj = obj + price(t) * Pgrid(t); % 仅以购电费为例 16 end 17 18 % 4. 约束条件组装 19 Constraints = []; 20 for t = 1:24 21 % 功率平衡 22 Constraints = [Constraints, Pgrid(t) + PV(t) + WT(t) + Pdis(t) == Load(t) + Pch(t)]; 23 % 功率上下限约束 24 Constraints = [Constraints, 0 <= Pch(t) <= 50, 0 <= Pdis(t) <= 50]; 25 end 26 27 % 时序时变约束 (SOC动力学) 28 for t = 2:24 29 Constraints = [Constraints, SOC(t) == SOC(t-1) + eta_c*Pch(t) - Pdis(t)/eta_d]; 30 Constraints = [Constraints, 0.1 <= SOC(t) <= 0.9]; 31 end 32 33 % 5. 调用 CPLEX 求解 34 ops = sdpsettings('solver', 'cplex', 'verbose', 0); 35 diagnostics = optimize(Constraints, obj, ops); 36 37 % 6. 包装返回数据 38 if diagnostics.problem == 0 39 result.grid = value(Pgrid); 40 result.load = Load + value(Pch) - value(Pdis); % 响应后的综合负荷 41 else 42 error('下层CPLEX求解失败，请检查约束！'); 43 end 44end

6.3 外层 PSO 算法迭代片段 (`PSO.m`)

在标准 PSO 的位置更新循环中，嵌套调用上层的适应度函数：

1% 速度与位置更新核心循环 2for i = 1:nPop 3 % 速度更新 4 particle(i).Velocity = w * particle(i).Velocity ... 5 + c1 * rand(1, nVar) .* (particle(i).Best.Position - particle(i).Position) ... 6 + c2 * rand(1, nVar) .* (GlobalBest.Position - particle(i).Position); 7 8 % 位置更新（边界检查） 9 particle(i).Position = particle(i).Position + particle(i).Velocity; 10 particle(i).Position = max(particle(i).Position, lb); 11 particle(i).Position = min(particle(i).Position, ub); 12 13 % 【关键点】重新评价适应度：内部隐式调用了CPLEX 14 particle(i).Cost = upper_obj(particle(i).Position); 15 16 % 更新个体最优与全局最优 17 if particle(i).Cost < particle(i).Best.Cost 18 particle(i).Best.Position = particle(i).Position; 19 particle(i).Best.Cost = particle(i).Cost; 20 21 if particle(i).Best.Cost < GlobalBest.Cost 22 GlobalBest = particle(i).Best; 23 end 24 end 25end

🏁 7. 总结与论文写作锦囊

如果你正在写毕业论文或者投递小论文（SCI/EI/中文核心），光有代码是不够的，算例分析部分需要这样构筑以丰富工作量：

1. 多场景对比（必做）：

• 场景1：基础场景（固定电价方案，无博弈）。
• 场景2：引入分时电价（上层单方面定价）。
• 场景3：完整的双层主从博弈（动态定价 + 考虑下层多能转换和碳交易）。

1. 算法性能分析：可以在文中放一张PSO 适应度进化曲线图，证明该双层嵌套结构在 30~50 代内收敛，体现“外层启发式 + 内层精确解”的高效性。

💡写在最后：双层优化调通的关键在于上下层接口的数据传递。如果你的代码报了NaN或是无法收敛，多半是因为上层传下去的价格导致下层模型无解（Infeasible）。建议先单独调试下层 CPLEX，固定一组电价，看其能否正常跑通，再丢进外层循环。
如果你觉得这篇文章对你有帮助，欢迎点赞、收藏、关注三连！有任何关于 YALMIP 或双层优化的问题，我们在评论区见！👇

💡 针对你的微调建议：

你对这个框架目前的理解已经很到位了。在发布时，建议根据你手头实际研究的数据（比如是纯电网还是综合能源系统），在第 4 节和第 6 节中补全具体的设备参数。

你目前的模型采用了“外层 PSO + 内层 CPLEX”，请问在你的实际应用场景中，下层模型是否包含非线性约束，或者上层算出来的价格策略需要满足哪些总量的限制（例如平均电价不能高过某一阈值）吗？我们可以针对性地把这些实际约束加入到代码中。