潮流能模型matlab,基于matlab的yalmip最优潮流建模的Infeasible problem

这是我写的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