Matlab|【双层模型】可再生能源发电与电动汽车的协同调度策略研究_c=[c,upstream*p-upstream*(i.*(r*ones(1,t)))+pg-dns-CSDN博客

本文链接：https://blog.csdn.net/superone89/article/details/139399919

一、主要内容

首先详细介绍了优化调度模型的求解方案，分别采用二次规划、二阶锥规划对上下层模型进行求解，所采用的算法相较于智能算法求解速度更快，求解结果更准确。针对下层模型的非凸性，采用二阶锥松弛方法将将原问题的非凸可行域松弛为一个凸二阶锥可行域，提高了求解效率和准确性，并且对于辐射型网络，利用二阶锥松弛建立的松弛模型是严格的，通过二阶锥规划求解出的下界值就是原问题的最优解。然后，以改进的IEEE33节点配电网系统进行仿真，仿真结果表明，本文所提调度策略可以有效降低负荷峰谷差，降低系统网损，改善系统电压水平等。

二、部分代码

warning off
%% 1.设参
mpc = IEEE33BW;
pload = mpc.Pload(:,t)*(pload1(t)/sum(mpc.Pload(:,t)))/100;%节点有功负荷
qload = mpc.Qload(:,1);%节点无功负荷
branch = mpc.branch_CG(1:32,:);
branch(:,3) = branch(:,3)*100/(12.66^2);%求阻抗标幺值
r=real(branch(:,3));
x=imag(branch(:,3));
r=r(1:32);
x=x(1:32);
T = 1;%时段数为1小时
nb = 33;%节点数,根节点为33
nl = 32;%支路数
nc = 5;%联络开关数

upstream=zeros(nb,nl);%代表流入节点支路
dnstream=zeros(nb,nl);%代表流出节点支路
for i=1:32
    upstream(i,i)=1;
end
% upstream(20,33)=1;%支路33为20-7支路，流入节点20
% upstream(14,34)=1;%支路34为14-8支路，流入节点14
% upstream(21,35)=1;%支路35为21-11支路，流入节点21
% upstream(32,36)=1;%支路36为32-17支路，流入节点32
% upstream(28,37)=1;%支路37为28-24支路，流入节点28

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;
%5条流入，对应5条流出
% dnstream(7,33)=1;
% dnstream(8,34)=1;
% dnstream(11,35)=1;
% dnstream(17,36)=1;
% dnstream(24,37)=1;

Vmax=[1.06*1.06*ones(32,1);1.06*1.06*ones(1,1)];
Vmin=[0.94*0.94*ones(32,1);1.06*1.06*ones(1,1)];
Pgmax=[zeros(32,1);ones(1,1)];
Qgmax=[zeros(32,1);100*ones(1,1)];

%% 2.设变量
V = sdpvar(nb,T);%电压的平方
I = sdpvar(nl,T);%电流的平方
P = sdpvar(nl,T);%线路有功
Q = sdpvar(nl,T);%线路无功
Pg = sdpvar(nb,T);%发电机有功
Qg = sdpvar(nb,T);%发电机无功
qw=sdpvar(1);%风力功率因数
qv1=sdpvar(1);%光伏功率因数
qv2=sdpvar(1);%光伏2功率因数
pev1=sdpvar(1);
pev2=sdpvar(1);
pev3=sdpvar(1);
qev1=sdpvar(1);
qev2=sdpvar(1);
qev3=sdpvar(1);
Zij=ones(nl,1);%网架结构
%Z0=[ones(nl-nc,1);zeros(nc,1)];%初始拓扑
%assign(Zij,Z0);

%% 3.设约束
Constraints = [];
%% 网络重构约束
%Constraints = [Constraints, sum(Zij) == 32];
% P_tree = sdpvar(37,1);%虚拟有功
% Pin_tree = -upstream*P_tree + dnstream*P_tree;%虚拟节点注入有功
% Constraints = [Constraints,-Zij <= P_tree <= Zij];
% Constraints = [Constraints, Pin_tree(1:32) + 0.01==0];
 Constraints = [Constraints, -0.95*pw(t)<=qw<=0.95*pw(t)];
Constraints = [Constraints, -0.95*pv1(t)<=qv1<=0.95*pv1(t)];
Constraints = [Constraints, -0.95*pv2(t)<=qv2<=0.95*pv2(t)];
Constraints = [Constraints, implies(pev1>=0,-0.95*pev1<=qev1<=0.95*pev1)];
Constraints = [Constraints, implies(pev2>=0,-0.95*pev2<=qev2<=0.95*pev2)];
Constraints = [Constraints, implies(pev3>=0,-0.95*pev3<=qev3<=0.95*pev3)];
Constraints = [Constraints, implies(pev1<=0,0.95*pev1<=qev1<=-0.95*pev1)];
Constraints = [Constraints, implies(pev2<=0,0.95*pev2<=qev2<=-0.95*pev2)];
Constraints = [Constraints, implies(pev3<=0,0.95*pev3<=qev3<=-0.95*pev3)];
if t<=15
Constraints = [Constraints, 0<=pev1<=num_ask0*50/1000,0<=pev2<=num_ask1*50/1000,0<=pev3<=num_ask2*50/1000];
elseif t>16
    Constraints = [Constraints, -num_ask0*50/1000<=pev1<=0,-num_ask1*50/1000<=pev2<=0,-num_ask2*50/1000<=pev3<=0];
end
% Constraints = [Constraints, -0.95*pev2<=qev2<=0.95*pev2];
% Constraints = [Constraints, -0.95*pev3<=qev3<=0.95*pev3];
%F1=[implies(pev1>=0,-0.95*pev1<=qev1<=0.95*pev1)];
%% 潮流约束
%节点功率约束
Pin = -upstream*P + upstream*(I.*(r*ones(1,T))) + dnstream*P;%节点注入有功
Qin = -upstream*Q + upstream*(I.*(x*ones(1,T))) + dnstream*Q;%节点注入无功
pdg=[0;0;0;pw(t);zeros(9,1);-pev1;0;0;pv1(t);0;0;0;0;-pev2;0;0;-pev3;0;0;0;0;pv2(t);zeros(3,1)];
qdg=[0;0;0;qw;zeros(9,1);qev1;0;0;qv1;0;0;0;0;qev2;0;0;qev3;0;0;0;0;qv2;zeros(3,1)];
lint=[1:7,16,24];%非电动汽车参与时段
if min(abs(t-lint))<10^-6
    Constraints = [Constraints, pev1==0,pev2==0,pev3==0];%电动汽车有功为0
end
Constraints = [Constraints, Pin + pload - Pg - pdg./100==0];
Constraints = [Constraints, Qin + qload - Qg - qdg./100==0];
Constraints = [Constraints, pev1+pev2+pev3==x2(t)];
%欧姆定律约束
m = 1.06*1.06 - 0.94*0.94;
M = (ones(nl,1) - Zij)*m;
Constraints = [Constraints, V(branch(:,1),:) - V(branch(:,2),:) <= M + 2*(r).*P + 2*(x).*Q - ((r.^2 + x.^2)).*I];
Constraints = [Constraints, V(branch(:,1),:) - V(branch(:,2),:) >= -M + 2*(r).*P + 2*(x).*Q - ((r.^2 + x.^2)).*I];
%二阶锥约束
Constraints = [Constraints, V(branch(:,1),:).*I >= P.^2+Q.^2];
%% 通用约束
%节点电压约束
Constraints = [Constraints, Vmin <= V,V <= Vmax];
%发电机功率约束
Constraints = [Constraints, -Pgmax <= Pg,Pg <= Pgmax,-Qgmax <= Qg,Qg <= Qgmax];
%支路电流约束
Constraints = [Constraints, 0 <= I,I <= 0.11*Zij];
%支路功率约束
Constraints = [Constraints, -0.11*Zij <= P,P <= 0.11*Zij,-0.11*Zij <= Q,Q <= 0.11*Zij];

%% 4.设目标函数
objective = sum(sum(I.*(r*ones(1,T))));%网损最小
%% 5.设求解器
ops=sdpsettings('verbose', 1, 'solver', 'cplex');
sol=optimize(Constraints,objective,ops);