MATLAB代码:考虑电动汽车时空调度的配电网潮流优化 参考文档:He L , Yang J , Yan J , et al. A bi-layer optimization based temporal and spatial scheduling for large-scale electric vehicles[J]. Applied Energy, 2016, 168(apr.15):179-192. DOI:10.1016/j.apenergy.2016.01.089 主要内容:机组组合采用原文相同的线性化方法最优潮流采用二阶锥松弛替代原文算法

前阵子帮小区物业调配电系统，物业大叔愁得头发都掉了几根——夏天晚上七八点，大家下班回家开空调、充电动车，配电箱直接过载跳闸，连路灯都暗了半截。其实这就是电动汽车并网调度最常见的痛点：用户想省钱找谷时段充电，但扎堆充电又会把电网搞崩，得有个法子平衡两边的需求。

MATLAB代码:考虑电动汽车时空调度的配电网潮流优化 参考文档:He L , Yang J , Yan J , et al. A bi-layer optimization based temporal and spatial scheduling for large-scale electric vehicles[J]. Applied Energy, 2016, 168(apr.15):179-192. DOI:10.1016/j.apenergy.2016.01.089 主要内容:机组组合采用原文相同的线性化方法最优潮流采用二阶锥松弛替代原文算法

刚好翻到2016年那篇发在Applied Energy上的经典论文，双层调度的思路挺绝，但原文用的下层迭代算法算起来太磨人，我这次换了二阶锥松弛来做最优潮流，再配上和原文一致的机组组合线性化套路，还写了个能直接跑的Matlab代码，今天就把踩过的坑和干货都唠出来。

先唠明白为啥要这么改

原文的双层框架其实很好懂：上层管配电网的机组启停和出力计划，下层管每台EV车主的充电需求，两边互相妥协找最优解。但原文下层用的非线性潮流求解不仅慢，还容易卡在局部最优，换成二阶锥松弛之后，就能把非线性的潮流方程转换成凸优化问题，求解器几秒就能出全局最优解，而且机组组合的线性化套路和原文完全对齐，不用改太多逻辑。

直接上代码，边写边拆

咱们用最简单的3节点配电网当测试案例，节点1是主网联入的Slack节点，节点2和3各装了5台充电桩，调度周期是24小时，每小时一个时段。首先初始化参数和工具箱：

%% 考虑EV调度的配电网潮流优化（二阶锥松弛版）
% 对齐Applied Energy 2016论文的机组组合逻辑，替换下层潮流为二阶锥松弛
clear;clc;close all
% 装了YALMIP和Gurobi的同学直接跑，没装的换intlinprog也能跑，就是慢

第一步：初始化系统参数

先把配电网、发电机、EV的基础参数都写清楚，不用太复杂，能跑通就行：

%% 1. 系统参数初始化
bus_num = 3;          % 总节点数
gen_num = 1;          % 配网内传统发电机组数量
ev_num_per_bus =5;    % 每个充电桩节点的EV数量
T = 24;               % 调度时段：24小时

% 发电机组参数
P_gen_max = 10;       % 最大出力MW
P_gen_min = 2;        % 最小出力MW
R_up = 3; R_down =3;  % 上下爬坡速率MW/h
fuel_cost = 0.1;      % 发电燃料成本 元/MWh

% EV参数
ev_cap =50;           % 单台EV电池容量kWh
ev_charge_max =3;     % 单台EV最大充电功率kW
% 车主默认下班回家要把电池从20%充到90%，单台总充电需求35kWh
ev_total_demand = 35*ev_num_per_bus/1000; % 转成MW·h

% 基础负荷（不含EV充电）
base_load = [5;3;4];  % 节点1是Slack，节点2、3的基础负荷
% 配网支路参数：支路1-2，1-3，阻抗r+jx
branch = [0.01 0.1; 0.02 0.15];
bus_from = [1;1]; bus_to = [2;3];
% 峰谷电价：白天8-22点电价高
price = ones(1,T)*0.5;
price(8:22) =1.2;

第二步：定义优化变量

这里用YALMIP的变量语法，比自己手写矩阵方便太多：

%% 2. 定义优化变量
% 发电机启停状态（二进制变量）和出力
u_gen = binvar(gen_num,T);
P_gen = sdpvar(gen_num,T);

% 二阶锥松弛专用变量：节点电压平方，避开极坐标的非线性sin/cos
V2 = sdpvar(bus_num,T);

% EV充电功率：节点2、3的总充电功率
P_ev_total = sdpvar(2,T); % 行对应节点2、3，列对应时段

% 支路有功/无功功率，用来写潮流约束
P_branch = sdpvar(size(branch,1),T);
Q_branch = sdpvar(size(branch,1),T);

第三步：添加约束，和原文对齐

首先是机组组合的线性化约束，和原文的逻辑完全一致：

%% 3. 添加约束
Constraints = [];
% -------------------------- 机组组合约束 --------------------------
for t =1:T
    % 发电机出力必须在启停状态对应的上下限之间
    Constraints = [Constraints, P_gen(:,t) >= P_gen_min*u_gen(:,t)];
    Constraints = [Constraints, P_gen(:,t) <= P_gen_max*u_gen(:,t)];
    % 爬坡约束：防止机组出力突变
    if t>1
        Constraints = [Constraints, P_gen(:,t)-P_gen(:,t-1) <= R_up*u_gen(:,t) + P_gen_max*(1-u_gen(:,t-1))];
        Constraints = [Constraints, P_gen(:,t-1)-P_gen(:,t) <= R_down*u_gen(:,t-1) + P_gen_max*(1-u_gen(:,t))];
    end
end
% 初始状态默认机组开机，避免凌晨突然停机
Constraints = [Constraints, u_gen(:,1)==1];

% -------------------------- EV充电约束 --------------------------
for b=1:2 % 只处理节点2、3的EV
    for t=1:T
        % 单节点总充电功率不能超过上限
        Constraints = [Constraints, P_ev_total(b,t) <= ev_num_per_bus*ev_charge_max/1000];
        Constraints = [Constraints, P_ev_total(b,t) >=0];
        % 全天总充电量必须满足车主需求
        if t==T
            Constraints = [Constraints, sum(P_ev_total(b,:)) >= ev_total_demand];
        end
    end
end

% -------------------------- 配网潮流约束（二阶锥松弛版） --------------------------
% Slack节点电压固定为1.0p.u.，平方就是1
Constraints = [Constraints, V2(1,:)==1];
% 节点电压必须在合格范围内：0.95~1.05p.u.
for t=1:T
    Constraints = [Constraints, 0.95^2 <= V2(:,t) <=1.05^2];
end

% 支路功率的二阶锥约束，还有节点有功平衡
for k=1:size(branch,1)
    % 计算支路导纳
    r = branch(k,1);x=branch(k,2);
    g = r/(r^2+x^2); b = -x/(r^2+x^2);
    for t=1:T
        Vi = sqrt(V2(bus_from(k),t));
        Vj = sqrt(V2(bus_to(k),t));
        % 二阶锥松弛核心：支路功率模长不超过电压平方的组合
        Constraints = [Constraints, norm([P_branch(k,t);Q_branch(k,t)]) <= sqrt(V2(bus_from(k),t)*V2(bus_to(k),t))*sqrt(g^2 +b^2)];
        
        % 节点有功平衡：总注入功率 = 基础负荷 + EV负荷 - 支路流出功率
        for b_node=2:bus_num
            if b_node ==1
                net_power = P_gen(1,t) - base_load(b_node) - P_ev_total(b_node-1,t);
            else
                net_power = -base_load(b_node) - P_ev_total(b_node-1,t);
            end
            % 统计所有连接到该节点的支路功率
            branch_flow =0;
            if bus_from(k)==b_node
                branch_flow = branch_flow + P_branch(k,t);
            end
            if bus_to(k)==b_node
                branch_flow = branch_flow - P_branch(k,t);
            end
            Constraints = [Constraints, net_power == branch_flow];
        end
    end
end

第四步：定义目标函数并求解

目标函数和原文一致，最小化总运行成本：包括发电燃料成本和EV充电的峰谷电价成本

%% 4. 目标函数与求解
% 总发电成本
gen_cost = fuel_cost * sum(sum(P_gen));
% 总EV充电成本：MW·h * 元/MWh = 元
ev_cost = sum(sum(P_ev_total.*repmat(price,2,1)));
Objective = gen_cost + ev_cost;

% 调用求解器，用Gurobi最快，没装的换'solver','intlinprog'
ops = sdpsettings('solver','gurobi','verbose',1);
result = optimize(Constraints,Objective,ops);

第五步：结果展示

跑通之后直接画图看结果：

%% 5. 结果展示
if result.problem ==0
    disp('恭喜！求解成功啦～');
    P_gen_res = value(P_gen);
    V2_res = value(V2);
    P_ev_res = value(P_ev_total);
    
    % 画发电机出力曲线
    figure('Name','机组出力');
    plot(1:T,P_gen_res(1,:),'linewidth',2);
    xlabel('时段(h)');ylabel('出力(MW)');title('传统机组出力');grid on;
    
    % 画EV充电曲线
    figure('Name','EV充电调度');
    plot(1:T,P_ev_res(1,:),'r',1:T,P_ev_res(2,:),'b','linewidth',2);
    xlabel('时段(h)');ylabel('总充电功率(MW)');legend('节点2','节点3');grid on;
    
    % 画节点电压
    figure('Name','节点电压');
    plot(1:T,sqrt(V2_res(2,:)),'r',1:T,sqrt(V2_res(3,:)),'b',1:T,ones(1,T),'k--','linewidth',1.5);
    xlabel('时段(h)');ylabel('电压(p.u.)');legend('节点2','节点3','标准电压');grid on;
else
    disp('坏啦，求解失败，错误码：');disp(result.problem);
end

踩过的坑唠两句

一开始我直接用非线性潮流写，结果求解器跑了十分钟都没出结果，换成二阶锥松弛之后秒出。不过要注意几个坑：

1. 二阶锥松弛对于径向配电网是精确的，不用担心里解不出来可行解，但如果是网状配网可能需要额外的约束；
2. 一开始忘了加电压上下限，结果算出来电压跑到1.5p.u.，明显不符合电网标准；
3. EV的充电约束一开始只加了单台功率上限，忘了总充电量要满足车主的需求，结果算出来EV根本没充满电。

最后补两句

这个代码只是个入门demo，要是真的做大规模配网调度，还得加上EV的到达/离开时间、储能装置、不同充电桩的功率限制这些细节，但作为理解双层调度和二阶锥松弛的例子已经够用了。要是把代码改成和原文完全一致的双层迭代结构，也只需要把下层的EV调度单独抽出来循环求解就行，感兴趣的同学可以自己改改~