1 主要内容

该程序复现博士文章《互动环境下分布式电源与电动汽车充电站的优化配置方法研究》第四章《考虑充电负荷空间可调度特性的分布式电源与电动汽车充电站联合配置方法》，本章构建了考虑充电负荷空间可调度特性的分布式电源与电动汽车充电站联合配置模型，以年化社会总成本最小为目标确定配电系统中分布式电源与电动汽车充电站的最佳安装位置和安装容量。该模型以线性化的Distflow潮流方程表征系统中各个状态量间的关系，并应用二阶锥模型进行求解。最后，本章将IEEE-33节点配电系统与江苏地区一个实际的地理区域相耦合，并以此为算例系统对光伏电站、微型燃气轮机、电动汽车充电站的最优配置方案进行了求解，验证了所构建模型的有效性。

目标函数

该部分目标函数考虑比较全面，因此也是复现的难点之一。

电动汽车负荷建模

算例系统图

程序亮点

本程序将33节点系统和实际地理图有机结合，充分考虑了电动车充电备选节点和调度具体的影响，最终规划得到分布式电源、充电桩等的配置结果，很有参考价值！

2 部分代码

Jz=inv(J);%求解雅可比矩阵逆矩阵
B4j=Jz(33:64,33:64);%求解得到B4j值,未考虑平衡节点
B4=[B4j,zeros(32,1);zeros(1,33)];%得到包含平衡节点的B4
%%
branch = mpc.branch;
branch(:,3) = branch(:,3)*100/(12.66^2);%求阻抗标幺值
r=real(branch(:,3));
x=imag(branch(:,3));
r=r(1:32);
x=x(1:32);
upstream=zeros(nb,nl);%代表流入节点支路
dnstream=zeros(nb,nl);%代表流出节点支路
for i=1:32
    upstream(i,i)=1;
end
for i=[1:16,18:20,22:23,25:31]
    dnstream(i,i+1)=1;
end
dnstream(1,18)=1;
dnstream(2,22)=1;
dnstream(5,25)=1;
dnstream(33,1)=1;

Vmax=[1.05*1.05*ones(32,1);1.05*1.05*ones(1,1)];
Vmin=[0.95*0.95*ones(32,1);1.05*1.05*ones(1,1)];
Pgmax=[zeros(32,1);100.*ones(1,1)];
Qgmax=[zeros(32,1);100.*ones(1,1)];
%计算负荷矩阵
Pload=busd(:,7);
Qload=busd(:,8);
P1=repmat(ploadz,33,1).*repmat(Pload,1,4*T)/100;%有功负荷矩阵
Q1=repmat(ploadz,33,1).*repmat(Qload,1,4*T)/100;%无功负荷矩阵
P2=zeros(33,4*T*N);
Q2=zeros(33,4*T*N);
for i=1:4*T%转化为4*T*N列矩阵
    for j=1:N
    P2(:,5*i+j-5)=P1(:,i);
    Q2(:,5*i+j-5)=Q1(:,i);
    end
end
%定义变量
V = sdpvar(nb,4*T*N);%电压的平方
I = sdpvar(nl,4*T*N);%电流的平方
P = sdpvar(nl,4*T*N);%线路有功
Q = sdpvar(nl,4*T*N);%线路无功
Pg = sdpvar(nb,4*T*N);%发电机有功
Qg = sdpvar(nb,4*T*N);%发电机无功
Npv = intvar(8,1);%光伏节点安装数量
%pv = sdpvar(nb,4*T*N);
qv = sdpvar(nb,4*T*N);
qv_h=sdpvar(8,4*T*N);
%sv = sdpvar(nb,4*T*N);
Ng = intvar(6,1);%燃气轮机安装数量
pg = intvar(nb,4*T*N);
Constraints = [];
%光伏处理
unitpv=10*1e-5;%单位光伏容量
pv_h=repmat(Npv,1,4*T*N).*repmat(dw_pv,8,4*T).*unitpv/10;
pv=[zeros(4,4*T*N);pv_h(1,:);zeros(5,4*T*N);pv_h(2,:);zeros(2,4*T*N);pv_h(3,:);zeros(1,4*T*N);pv_h(4,:);zeros(3,4*T*N);pv_h(5,:);zeros(2,4*T*N);pv_h(6,:);zeros(5,4*T*N);pv_h(7,:);zeros(1,4*T*N);pv_h(8,:);zeros(2,4*T*N)];
sv_h=unitpv.*repmat(Npv,1,4*T*N);
sv=[zeros(4,4*T*N);sv_h(1,:);zeros(5,4*T*N);sv_h(2,:);zeros(2,4*T*N);sv_h(3,:);zeros(1,4*T*N);sv_h(4,:);zeros(3,4*T*N);sv_h(5,:);zeros(2,4*T*N);sv_h(6,:);zeros(5,4*T*N);sv_h(7,:);zeros(1,4*T*N);sv_h(8,:);zeros(2,4*T*N)];

% Constraints = [Constraints,pv>=0,pv<=repmat(unitpv.*Npv,1,4*T*N)];
% Constraints = [Constraints,pv>=0,pv<=repmat(unitpv.*Npv,1,4*T*N)];
for i=1:33
    for t=1:T
        Constraints = [Constraints,cone([pv(i,t);qv(i,t)],sv(i,t))];
    end
end
qv=[zeros(4,4*T*N);qv_h(1,:);zeros(5,4*T*N);qv_h(2,:);zeros(2,4*T*N);qv_h(3,:);zeros(1,4*T*N);qv_h(4,:);zeros(3,4*T*N);qv_h(5,:);zeros(2,4*T*N);qv_h(6,:);zeros(5,4*T*N);qv_h(7,:);zeros(1,4*T*N);qv_h(8,:);zeros(2,4*T*N)];
Constraints = [Constraints,Npv>=0,Ng>=0,Npv<=100];
%微燃机处理