这是我写的ieee33算例静态配电网潮流模型,使用了cplex,sedumi,mosek求解器均为Infeasible problem,请大佬为小白答疑解惑,应该如何解决呢?
clc
clear all
%% 1.设参
Gen_Pmax = 10; %发电机发出有功上限(MW)
Gen_Pmin = 0; %发电机发出有功下限
Gen_Qmax = 10; %发电机发出无功上限(MW)
Gen_Qmin = - Gen_Qmax; %发电机发出无功下限
bus_Vmax = 1.1^2; %节点电压上限
bus_Vmin = 0.95^2; %节点电压下限
bus_num = 33; %母线数
branch_num = 32; %支路数
Gen_num = 1; %发电机数
PLoadBus = [0 0.10 0.09 0.12 0.06 0.06 0.20 0.20 0.06 0.06 0.04 0.06 0.06 0.12 0.06 0.06 0.06 0.09 0.09 0.09 0.09 0.09 0.09 0.42 0.42 0.06 0.06 0.06 0.12 0.20 0.15 0.21 0.06];
QLoadBus = [0 0.06 0.04 0.08 0.03 0.02 0.10 0.10 0.02 0.02 0.03 0.035 0.035 0.08 0.01 0.02 0.02 0.04 0.04 0.04 0.04 0.04 0.05 0.20 0.20 0.025 0.025 0.02 0.07 0.60 0.07 0.10 0.04];
Branch = [1 2 0.0922 0.0470; 2 3 0.4930 0.2511;3 4 0.3660 0.1864; 4 5 0.3811 0.1941; 5 6 0.8190 0.7070;6 7 0.1872 0.6188;
7 8 0.7114 0.2351; 8 9 1.0300 0.7400; 9 10 1.0440 0.7400; 10 11 0.1966 0.0650; 11 12 0.3744 0.1238; 12 13 1.4680 1.1550;
13 14 0.5416 0.7129; 14 15 0.5910 0.5260; 15 16 0.7463 0.5450; 16 17 1.2890 1.7210; 17 18 0.3720 0.5740; 2 19 0.1640 0.1565;
19 20 1.5042 1.3554; 20 21 0.4095 0.4784; 21 22 0.7089 0.9373; 3 23 0.4512 0.3083; 23 24 0.8980 0.7091; 24 25 0.8960 0.7011;
6 26 0.2030 0.1034; 26 27 0.2842 0.1447; 27 28 1.0590 0.9337; 28 29 0.8042 0.7006; 29 30 0.5075 0.2585; 30 31 0.9744 0.9630;
31 32 0.3105 0.3619; 32 33 0.3410 0.5362];
%% 2.有关节点功率方程的矩阵
% 节点-发电机关联矩阵
GenBus = 1;
Gconnect = zeros(bus_num,Gen_num);
for i = 1:Gen_num
Gconnect(GenBus(i),i)=1;
end
%节点-支路关联矩阵(定义流出为正,流入为负)
Bconnect = zeros(bus_num, branch_num);
for i = 1 : branch_num
Bconnect(Branch(i, 1), i) = 1;
Bconnect(Branch(i, 2), i) = -1;
end
%节点-阻抗关联矩阵(由节点-支路关联矩阵改变而来)
Rconnect = zeros(bus_num, branch_num);
Yconnect = zeros(bus_num, branch_num);
for i = 1 : branch_num
Rconnect(Branch(i, 2), i) = Branch(i, 3);
Yconnect(Branch(i, 2), i) = Branch(i, 4);
end
%% 3.设决策变量P,Q,p,q,V,I(此V,I为相角松弛后的v,i,即V=v^2,I = i^2)
Pij = sdpvar(branch_num, 1, 'full');%x1-32
Qij = sdpvar(branch_num, 1, 'full'); %x33-64
P_Gen = sdpvar(Gen_num,1, 'full');%x65
Q_Gen = sdpvar(Gen_num,1, 'full');%x66
Vij = sdpvar(bus_num, 1, 'full');%x67-98
Iij = sdpvar(branch_num, 1, 'full');%x99-130
%% 4.设约束
Constraints = [];
Constraints = [Constraints, Gconnect*P_Gen == Bconnect*Pij + Rconnect*Iij + PLoadBus'];
Constraints = [Constraints, Gconnect*Q_Gen == Bconnect*Qij + Yconnect*Iij + QLoadBus'];
% Ohm's law约束
for i = 1:branch_num
Constraints = [Constraints, Vij(Branch(i, 2)) == Vij(Branch(i, 1)) - 2*(Pij(i)*Branch(i, 3) + Qij(i)*Branch(i, 4)) + Iij(i)*(Branch(i, 3)^2 + Branch(i, 4)^2)];
end
%使用norm
% m = [2*Pij;2*Qij;Iij-Vij(Branch(:, 1))];
% Constraints = [Constraints, norm(m) <= Iij+Vij(Branch(:, 1))];
%使用cone
for i = 1:branch_num
Constraints = [Constraints, cone([2*Pij(i);2*Qij(i);Iij(i)-Vij(Branch(i, 1))],Iij(i)+Vij(Branch(i, 1)))];
end
% 节点电压约束
Constraints = [Constraints, Vij(1) == 1];
for i = 1:bus_num
Constraints = [Constraints, bus_Vmin <= Vij(i) <= bus_Vmax];
end
%发电机功率约束
for i = 1:Gen_num
Constraints = [Constraints, Gen_Pmin <= P_Gen(i) <= Gen_Pmax];
Constraints = [Constraints, Gen_Qmin <= Q_Gen(i) <= Gen_Qmax];
end
%% 5.设目标函数
objective = 0;
for i = 1 : branch_num
objective = objective + Iij(i)*Branch(i, 3);
end
%% 6.设置求解器并求解
ops = sdpsettings('verbose', 1, 'solver', 'quadprog');
sol = optimize(Constraints, objective, ops);
Vij = value(Vij)
Pij = value(Pij)
%% 7.输出AMPL模型并设置行编号
saveampl(Constraints,objective,'mymodel');
%% 8.分析错误标志
if sol.problem == 0
disp('succcessful solved');
else
disp('error');
yalmiperror(sol.problem)
end