以上分析了综合负荷和部分需求响应,以下程序仅有综合负荷。
%% 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