【完整版代码】含分布式电源的配电网日前两阶段优化调度模型

目录

1 主要内容

1.1 摘要

1.2 具体模型

1.3 网架结构

2 部分程序

3 程序结果

4 程序链接

---

1 主要内容

该程序复现了《含分布式电源的配电网日前两阶段优化调度模型》-孟晓丽完整模型，包括第 1 阶段为 DG 优化调度阶段和第 2 阶段为无功优化阶段，之前发布的复现第二阶段程序链接如下：

《含分布式电源的配电网日前两阶段优化调度模型》第二阶段模型。

1.1 摘要

在电力市场环境下，供电公司通过对接入配电网的分布式电源(distributed generation，DG)的优化调度，能够有效地降低其运行成本，规避市场竞争环境下的风险。提出了一种电力市场环境下供电公司日前优化调度的 2 阶段模型：第 1 阶段为 DG 优化调度阶段，根据市场电价、DG运行成本、签订可中断负荷(interruptable load，IL)合同的价格来确定 DG 的机组组合、从大电网的购电量及 IL 削减量；第 2 阶段为无功优化阶段，在第 1 阶段的基础上，考虑 DG 的无功出力特性，通过优化 DG 和无功补偿装置的出力调节电压使其在规定的范围内且配电网的网损最小。通过基于修改的 IEEE 33 节点系统的仿真计算，表明所提出的日前 2 阶段优化调度模型能够有效降低供电公司的运行成本。

1.2 具体模型

• 第一阶段目标函数：

• 第一阶段约束条件：

有功平衡约束

系统备用约束

DG出力、爬坡、最小启停时间约束

大电网购电约束

IL合同约束

【思考】在复现第一阶段模型过程中，遇到几点疑问以及处理方式如下：

1.算例写明标准IEEE33节点配电系统总负荷为5084.26+j2547.32kVA

一般更常用的是总有功为3715kW的IEEE33节点数据文件，因此对负荷数据部分按照比例进行了调整。

2.系统备用约束增加了Rt作为新的变量，文章并没有给出Rt的其他约束形式，相当于这个约束失去了意义。

3.该模型忽略了一个关键的问题，本来为配电网模型，但是第一阶段完全是按照微网模型表达的，没有计及二阶锥约束、电压电流约束等配电网特征的约束，为了尽量靠近33节点模型，程序采用二阶锥约束形式的配电网模型。

4.一个小问题，图3有功功率单位标注出错。

5.结果图5和图6两部分功率加起来已经严重超过系统总负荷8MW(见图3），再叠加分布式电源出力，总供应功率远大于负荷功率，结果部分准确性存疑。

• 第二阶段目标函数如下：即网损和节点电压越限罚函数之和。

• 第二阶段约束条件如下：

潮流约束

节点电压约束

DG 无功补偿约束

无功补偿器容量约束

对于第 2 阶段，含分布式电源的配电网无功优化问题目标复杂、约束多，是一个多维的非线性优化问题，本文采用粒子群算法(particle swarm optimization，PSO)来求解。算法主要由潮流计算和优化两部分组成。其中潮流计算为优化计算提供状态变量的值和网损；优化部分主要完成粒子群算法中每代粒子的更新。

1.3 网架结构

本程序研究的33节点网架结构如下：

2 部分程序

x_Iij_square=sdpvar(32,N,'full');%电流平方  
x_ui_square=sdpvar(33,N,'full');%电压平方
x_pij=sdpvar(32,N,'full');%线路有功功率      
x_qij=sdpvar(32,N,'full');%线路无功功率
x_Pdg=sdpvar(3,N,'full');%dg功率
u_Pdg=binvar(3,N,'full');%dg出力状态
x_Pil=sdpvar(2,N,'full');%切负荷功率
u_Pil=binvar(2,N,'full');%切负荷状态
x0_ug=binvar(3,N,'full');%机组启动
x0_vg=binvar(3,N,'full');%机组停止
x_R=sdpvar(1,N,'full');%备用负荷
%% 定义约束条件
Constraints=[];
P_dg=sdpvar(33,24,'full');%将dg出力情况放入对应的节点
for i=1:33
    if i==7
    Constraints=[Constraints,P_dg(i,:)==x_Pdg(1,:)];
    elseif i==24
        Constraints=[Constraints,P_dg(i,:)==x_Pdg(2,:)];
    elseif i==25
        Constraints=[Constraints,P_dg(i,:)==x_Pdg(3,:)];
    else
        Constraints=[Constraints,P_dg(i,:)==0];
    end
end
P_il=sdpvar(33,24,'full');%将IL合同约束放入对应的节点
for i=1:33
    if i==8
    Constraints=[Constraints,P_il(i,:)==x_Pil(1,:)];
    elseif i==25
        Constraints=[Constraints,P_il(i,:)==x_Pil(2,:)];
    else
        Constraints=[Constraints,P_il(i,:)==0];
    end
end
Constraints=[Constraints,x_ui_square(1,:)==12.66^2];  %平衡节点每小时电压平方
Constraints=[Constraints,x_Iij_square>=0]; 
% Constraints=[Constraints,x_pij(1,:)>=0];
%% 支路欧姆定律
for r=1:32
    Constraints=[Constraints,x_ui_square(Branch(r,2),:)-...
        x_ui_square(Branch(r,3),:)+(r_ij(r)^2+x_ij(r)^2)*x_Iij_square(r,:)-...
        2*(r_ij(r)*x_pij(r,:)+x_ij(r)*x_qij(r,:))==0];
end
%dg爬坡约束
for t=1:23
    for i=1:3
    Constraints=[Constraints,x_Pdg(i,t+1)-x_Pdg(i,t)<=vup];
    Constraints=[Constraints,x_Pdg(i,t)-x_Pdg(i,t+1)<=vdn];
    end
end
%IL合同约束
Constraints=[Constraints,0<=x_Pil<=Pilmax.*u_Pil];
%dg最小启停时间约束
for i=1:3
Constraints=[Constraints,consequtiveON([onoffhistory(i,:) u_Pdg(i,:)],Tmup)];
Constraints=[Constraints,consequtiveON(1-[onoffhistory(i,:) u_Pdg(i,:)],Tmd)];
end
for t=2:N%公式7

3 程序结果