主动配电网SOCP_OPF学习笔记（5）综合负荷及其程序

       以上分析了综合负荷和部分需求响应，以下程序仅有综合负荷。

%% PG69最优潮流SOCP-OPF,SB=10MVA,VB=12.66kV,添加综合负荷
%恒功率 1-15，42-54
%恒电流 16-27，36-41，55-58
%恒阻抗 28-35，59-69
%%
clear
clc
tic
warning off
%% 1.获取数据
mpc = PG69BW;
Line = mpc.Line;
pload = 0.0001*[0;real(Line(:,4))];%转成MW并求标幺值
qload = 0.0001*[0;imag(Line(:,4))];
%负荷曲线形状系数*节点负荷=全天动态数据
Day = mpc.Day;
pload = pload*Day;%获取动态数据
qload = qload*Day;
SB = 10;
VB = 12.66;
ZB = VB^2/SB;
Line(:,3) = Line(:,3)/ZB;
r = real(Line(:,3));%阻抗标幺值
x = imag(Line(:,3));
%% 2.设参
nb = 69;
nl = 68;
T = 24;

%总时段潮流值
Ploss = zeros(T,1);%总损耗
V_total = zeros(nb,T);%节点电压
I_total = zeros(nl,T);%支路电流
P_total = zeros(nl,T);%支路有功
Q_total = zeros(nl,T);%支路无功
% Theta = zeros(nb,T);%节点相角

%% 3.潮流计算
parfor t= 1:T 
    %% ①设单时段潮流值
    V = sdpvar(nb,1);
    I = sdpvar(nl,1);
    P = sdpvar(nl,1);
    Q = sdpvar(nl,1);
    Pload = sdpvar(nb,1);%综合负荷有功
    Qload = sdpvar(nb,1);
    %% 综合负荷
    AA = [ones(15,1);zeros(26,1);ones(13,1);zeros(15,1)];%恒阻抗比例
    BB = [zeros(15,1);ones(12,1);zeros(8,1);ones(6,1);zeros(13,1);ones(4,1);zeros(11,1)];%恒电流比例
    CC = [zeros(27,1);ones(8,1);zeros(23,1);ones(11,1)];%恒功率比例
    Pload(:,t) = (AA+BB/2).*V.*pload(:,t)+(BB/2+CC).*pload(:,t);
    Qload(:,t) = (AA+BB/2).*V.*qload(:,t)+(BB/2+CC).*qload(:,t);
    %% ②设约束
    C=[];
    %节点功率约束
    for i = [1,4:6,9,12:25,27:33,35:37,39,41:52,54,56,58:67]
        C = [C, P(i) - I(i)*r(i) == Pload(i+1,t) + P(i+1)];
        C = [C, Q(i) - I(i)*x(i) == Qload(i+1,t) + Q(i+1)];
    end
    %多出线度节点
    C=[C, P(2) - I(2)*r(2) == Pload(3,t) + P(3) + P(58) + P(27)];
    C=[C, Q(2) - I(2)*x(2) == Qload(3,t) + Q(3) + Q(58) + Q(27)];
    C=[C, P(3) - I(3)*r(3) == Pload(4,t) + P(4) + P(35)];
    C=[C, Q(3) - I(3)*x(3) == Qload(4,t) + Q(4) + Q(35)];
    C=[C, P(7) - I(7)*r(7) == Pload(8,t) + P(8) + P(39)];
    C=[C, Q(7) - I(7)*x(7) == Qload(8,t) + Q(8) + Q(39)];
    C=[C, P(8) - I(8)*r(8) == Pload(9,t) + P(9) + P(41)];
    C=[C, Q(8) - I(8)*x(8) == Qload(9,t) + Q(9) + Q(41)];
    C=[C, P(10) - I(10)*r(10) == Pload(11,t) + P(11) + P(54)];
    C=[C, Q(10) - I(10)*x(10) == Qload(11,t) + Q(11) + Q(54)];
    C=[C, P(11) - I(11)*r(11) == Pload(12,t) + P(12) + P(56)];
    C=[C, Q(11) - I(11)*x(11) == Qload(12,t) + Q(12) + Q(56)];
    %辐射支路尾端节点
    for i = [26,34,38,40,53,55,57,68]
        C = [C, P(i) - I(i)*r(i) == Pload(i+1,t)];
        C = [C, Q(i) - I(i)*x(i) == Qload(i+1,t)];
    end
    %支路欧姆定律约束
    for i = 1:nl
        C = [C, V(i+1) == V(Line(i,1)) - 2*(P(i)*r(i) + Q(i)*x(i)) + (r(i)^2 + x(i)^2)*I(i)];
    end
    %SOC约束
    for i = 1:nl
        C = [C, V(Line(i,1))*I(i) >= P(i)^2 + Q(i)^2];
    end

    %节点电压约束
    C = [C, 0.94^2 <= V,V <= 1.06^2];
    V(1) = 1.06^2;
    %支路电流约束
    C = [C, 0 <= I,I <= 1];%未添加主动管理，支路电流不设限

    %% ③目标函数
     objective=sum(I.*r);
     
    %% ④求解
    ops = sdpsettings('solver','cplex','verbose',0);
    sol = optimize(C,objective,ops);
    
    %% ⑤保存单时间潮流
    Ploss(t) = double(objective);
    V_total(:,t) = sqrt(double(V));
    I_total(:,t) = sqrt(double(I));
    P_total(:,t) = double(P);
    Q_total(:,t) = double(Q);
    %% ⑥分析错误标志
    if sol.problem == 0
       disp('succcessful solved');
    else
       disp('error');
       yalmiperror(sol.problem)
     end
end
Ploss_total = 10*sum(Ploss);
clear C f i I Line ops P Pload Q Qload r result SB U UB x ZB Day t nb nl mpc Ploss VB T pload qload
toc