Yalmip使用教程(8)-常见报错及调试方法

本文主要是介绍Yalmip使用教程(8)-常见报错及调试方法，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

博客中所有内容均来源于自己学习过程中积累的经验以及对yalmip官方文档的翻译：https://yalmip.github.io/tutorials/

这篇博客将详细介绍使用yalmip工具箱编程过程中的常见错误和相应的解决办法。

1.optimize的输出参数

众所周知，optimize是yalmip用来求解优化问题的函数，其使用格式为：

sol = optimize(Constraints,Objective,options)

optimize函数的返回值sol是一个包含六个字段的结构体：

>> solsol = 包含以下字段的 struct:yalmipversion: '20210331'matlabversion: '9.1.0.441655 (R2016b)'yalmiptime: 0.0851solvertime: 0.0439info: 'Successfully solved (GUROBI-GUROBI)'problem: 0

其中，yalmipversion表示Yalmip工具箱的版本，matlabversion表示Matlab的版本，yalmiptime表示Yalmip的建模时间，solvertime表示求解器的求解时间，info表示返回的信息，problem为求解成功的标志，0表示求解成功，1表示求解失败。

其中最重要的参数就是problem和info，可以显示求解是否成功，以及可能遇到的问题。因此通常可以在optimize函数求解之后再加一部分代码来展示是否求解成功和求解失败的原因：

if sol.problem == 0disp('求解成功')
else disp('求解失败，失败原因为：')disp(sol.info)
end

以下的代码是一个能求解成功的例子：

x = sdpvar(2,1);
Constraints = [x <= 3 , x(1) >= 2 , x(2) >= 1];
Objective = [1 3]*x;
ops = sdpsettings;
sol = optimize(Constraints , Objective , ops);
if sol.problem == 0disp('求解成功')
else disp('求解失败，失败原因为：')disp(sol.info)
end

下面的例子是一个求解失败的代码：

x = sdpvar(2,1);
Constraints = [x <= 1 , x(1) >= 2 , x(2) >= 1];
Objective = [1 3]*x;
ops = sdpsettings;
sol = optimize(Constraints , Objective , ops);
if sol.problem == 0disp('求解成功')
else disp('求解失败，失败原因为：')disp(sol.info)
end

因为约束条件中x1≤1且x1≥2，所以导致优化问题无解。

在确保优化问题求解成功的情况下，可以采用value命令求出目标函数或者决策变量的取值，例如：

x = sdpvar(2,1);
Constraints = [x <= 3 , x(1) >= 2 , x(2) >= 1];
Objective = [1 3]*x;
ops = sdpsettings;
sol = optimize(Constraints , Objective , ops);
if sol.problem == 0disp('求解成功')x1 = value(x(1))x2 = value(x(2))Objective = value(Objective)
elsedisp('求解失败，失败原因为：')disp(sol.info)
end

结果为：

表示这个优化问题的最优解为x1=2，x2=1，目标函数最小值为5。

2.根据info信息判断判断解决方法

由第一节可知，yalmip工具箱中出现的大部分错误的原因都可以从optimize输出参数的info属性中查看。

例如，下面这个图是非常常见的错误：

主要原因是优化问题无解，optimize函数求解失败，后续又用到了优化问题中的一些变量，就会出现图中的报错。解决方法就是需要确保optimize函数可以求解成功。下面分别介绍一些常见的错误以及相应的解决方法。

2.1 solver not found

这个错误很简单，就是字面意思（没有安装求解器），也是私信问我的朋友中最常见的错误。

解决这个报错也很简单，没有求解器的去安装一个，或者换成已经安装过的求解器就行了。

和这个报错同类型的提示还有：

Solver license cannot be located

Solver license expired

Specified solver name not recognized

License problems in solver

上述提示均可采用相同的方式解决。

2.2 Solver not applicable

从字面意思看，这个提示就是说求解器不适用。可能这么说大家还是不太明白。举个例子，LINPROG是matlab自带的线性规划求解器，不具备求解二次规划的功能，如果我们用这个求解器来解二次规划问题，yalmip就会报错提示求解器不可用，如下面的代码：

clc
clearx1 = sdpvar(1);
x2 = sdpvar(1);
Constraints = [x1 >= 1 , x1 <= 2 , x2 <= 5 ,x2 >= 1];
Objective = x1^2 +9*x2;
ops = sdpsettings('solver' , 'LINPROG' , 'verbose' , 3 );
sol = optimize(Constraints , -Objective , ops);
sol.info

输出结果为：

ans ='Solver not applicable (linprog does not support nonconvex quadratic terms in objective)'

如果我们在求解优化问题的过程中出现上面的提示，就需要确定优化问题的类型，仔细检查你所使用的求解器是否支持该类型优化问题。换一个能求解该类型问题的求解器即可。

比如gurobi可以用来求解二次规划问题，我们可以把上面例子中的求解器改为gurobi：

clc
clearx1 = sdpvar(1);
x2 = sdpvar(1);
Constraints = [x1 >= 1 , x1 <= 2 , x2 <= 5 ,x2 >= 1];
Objective = x1^2 +9*x2;
ops = sdpsettings('solver' , 'gurobi' , 'verbose' , 3 );
sol = optimize(Constraints , -Objective , ops);
sol.info

这时候就可以求解成功了：

ans ='Successfully solved (GUROBI-NONCONVEX)'

2.3 Unbounded objective function

从字面意思看，这个提示就是指优化模型的目标函数是无界的。具体来说，就是目标函数中存在没有边界范围的变量，使得目标函数无法取得最大值（或最小值）。例如下面的优化问题：

很明显这是一个无界的优化问题，因为变量x没有下界，因此无法找到2x的最小值，我们用yalmip求解试试：

clc
clearsdpvar x
Constraints = [x <= 1];
Objective = 2*x;
ops = sdpsettings('solver' , 'cplex' , 'verbose' , 3 );
sol = optimize(Constraints , Objective , ops);
sol.info

结果如下：

ans ='Unbounded objective function (CPLEX-IBM)'

针对目标函数无界的情况，调试的方法也很简单。我们可以给定目标函数中涉及到的所有变量上下界，那么优化问题一般都可以从无解转为有解，如果存在某些变量，无论上下界取何值，该变量的取值都处于我们给定的边界上，那么就是因为这个变量没有给定范围，才导致目标函数无界。举个例子：

clc
clearsdpvar x y zConstraints = [-1 <= [x y] <= 1, z <= 0];
Objective = x^2 + y + z;
ops = sdpsettings('solver' , 'cplex' , 'verbose' , 3 );
sol = optimize(Constraints , Objective , ops);
sol.info

这个优化问题也是无界，为了找到具体是哪个变量导致目标函数无界，我们可以按下面的步骤操作：

1）用allvariables函数取出目标函数中所有涉及的变量（注意，使用allvariables函数要求yalmip版本高于2021，旧版工具箱里面不存在该函数，旧版yalmip可以手动选取目标函数涉及的所有变量）：

UsedInObjective = allvariables(Objective);

如果使用了较低版本的yalmip，可以把代码改成下面这样：

UsedInObjective = [x, y, z];

2）人为给定这些变量的上下限：

Constraints = [Constraints, -100 <= UsedInObjective <= 100];

3）求解修改后的优化问题，并查看所有变量的取值：

ops = sdpsettings('solver' , 'cplex' , 'verbose' , 3 );
sol = optimize(Constraints , Objective , ops);
value([x y z])

结果如下：

ans =0    -1  -100

4）修改这些变量的上下限，并重复步骤2和3：

Constraints = [Constraints, -1000 <= UsedInObjective <= 1000];
ops = sdpsettings('solver' , 'cplex' , 'verbose' , 3 );
sol = optimize(Constraints , Objective , ops);
value([x y z])

结果如下：

ans =0    -1  -1000

我们可以看到，无论我们给定多大的边界，变量z取值都位于边界，因此就是未定义变量z的范围或相关的约束，导致模型的目标函数无界，给定变量z的范围或相关的约束即可解决该问题。

2.4 Infeasible problem

这个提示的字面意思是“不可行的问题”，也就是优化问题是无解的。一般情况下，优化问题无解都是因为约束条件之间存在矛盾或者限制的太死，我们需要通过调试找到问题所在。对于大规模的优化问题，调试的过程是很痛苦的，我们在调试之前可以先做一些简单的检查。

2.4.1 一些简单检查

1）检查变量的定义是否存在问题，如0-1变量是否用binvar函数定义，连续变量是否用sdpvar函数定义，整数变量是否用intvar函数定义。

2）对于多维变量，检查是否添加了’full’属性。由于yalmip在定义N×N维变量时，如果不添加’full’属性，将默认变量是对称的，定义多维变量时很容易忽视这一点，导致模型不可行。

举个例子，我在之前的博客中提到了一个数独问题（Yalmip入门教程(3)-约束条件的定义-CSDN博客）：

求解该问题的代码如下：

clc
clearA0 = [
0 4 7 0 5 0 0 0 8;
6 0 5 0 3 0 2 0 1;
0 0 0 7 0 6 0 3 0;
0 0 6 0 7 0 0 2 4;
9 0 0 8 0 4 0 0 6;
4 5 0 0 1 0 9 0 0;
0 1 0 5 0 2 0 0 0;
2 0 8 0 4 0 5 0 3;
5 0 0 0 9 0 7 1 0;
];A = binvar(9,9,9,'full');
F = [sum(A,1) == 1, sum(A,2) == 1, sum(A,3) == 1];for i = 1:9for j = 1:9if A0(i,j)F = [F, A(i,j,A0(i,j)) == 1];endend
endfor m = 1:3for n = 1:3for k = 1:9s = sum(sum(A((m-1)*3+(1:3),(n-1)*3+(1:3),k)));             F = [F, s == 1];endend
enddiagnostics = optimize(F);Z = 0;
for i = 1:9Z = Z  + i*value(A(:,:,i));
end
Z

运行结果如下：

Z =3     4     7     2     5     1     6     9     86     8     5     4     3     9     2     7     11     2     9     7     8     6     4     3     58     3     6     9     7     5     1     2     49     7     1     8     2     4     3     5     64     5     2     6     1     3     9     8     77     1     3     5     6     2     8     4     92     9     8     1     4     7     5     6     35     6     4     3     9     8     7     1     2

如果我们将变量A定义时的’full’属性删除，再次运行的结果如下：

3）确定模型是否真的无解，而不是目标函数无界。有时候求解器会提示Either infeasible or unbounded，无法确定是优化问题无解还是目标函数无界。这时候我们可以将目标函数设为空，如果修改后可以求解成功，说明是目标函数无界，按照本博客2.3小节提到的方法进行调试即可。如果修改后还是提示模型不可行，那么就需要我们继续调试。

2.4.2 调试较大规模的不可行问题

简单来说，调试较大规模的不可行问题的思路就是依次向优化问题中添加约束条件，确定具体是哪组约束条件导致模型不可行或者哪些约束条件存在矛盾导致模型不可行。

我之前的博客给出了一个基于混合整数二阶锥(MISOCP)的配电网重构代码(基于混合整数二阶锥(MISOCP)的配电网重构(附matlab代码)-CSDN博客)，这是一个稍微有一丢丢复杂的优化问题，我悄咪咪地修改了其中一些参数，使得模型不可行，下面以修改后优化问题为例，讲一下调试不可行问题的思路。

首先给出修改后的代码：

%% 清除内存空间
clc
clear
close all
warning off
%% 系统参数
mpc = IEEE33;
nb=33;                              % 节点数
ns=1;                               % 电源节点数
nl=37;                              % 支路数
P_load=mpc.bus(:,3)/mpc.baseMVA;  % 有功负荷
Q_load=mpc.bus(:,4)/mpc.baseMVA;  % 无功负荷
r_ij=mpc.branch(:,3);              % 线路电阻
x_ij=mpc.branch(:,4);              % 线路电抗
M=1.06*1.06 - 0.94*0.94;
% 电源功率最大值
P_g_max=zeros(nb,1);
P_g_max(1)=10/mpc.baseMVA;
Q_g_max=zeros(nb,1);
Q_g_max(1)=10/mpc.baseMVA;
% 电源功率最小值
P_g_min=zeros(nb,1);
Q_g_min=[-10/mpc.baseMVA;zeros(nb-1,1)];
% 支路功率最大值
Sij_max=15/mpc.baseMVA;
% 最大电压
Umax=[1;1.1*1.1*ones(32,1)];
Umin=[1;0.95*0.95*ones(32,1)];
% 流入节点的支路
branch_to_node=zeros(nb,nl);
% 流出节点的支路
branch_from_node=zeros(nb,nl);
for k=1:nlbranch_to_node(mpc.branch(k,2),k)=1;branch_from_node(mpc.branch(k,1),k)=1;
end%% 优化变量
z_ij=binvar(nl,1);% 支路开断情况
U_i=sdpvar(nb,1);% 电压的平方
L_ij=sdpvar(nl,1);% 电流的平方
P_ij=sdpvar(nl,1);% 线路有功功率
Q_ij=sdpvar(nl,1);% 线路无功功率
P_g=sdpvar(nb,1);% 电源有功出力
Q_g=sdpvar(nb,1);% 电源无功功率%% 约束条件
Constraints = [];
%% 1.潮流约束
m_ij=(1-z_ij)*M; 
Constraints = [Constraints, P_g-P_load+branch_to_node*P_ij-branch_to_node*(L_ij.*r_ij)-branch_from_node*P_ij == 0];
Constraints = [Constraints, Q_g-Q_load+branch_to_node*Q_ij-branch_to_node*(L_ij.*x_ij)-branch_from_node*Q_ij == 0];
Constraints = [Constraints,U_i(mpc.branch(:,1))-U_i(mpc.branch(:,2))<= m_ij + 2*r_ij.*P_ij + 2*x_ij.*Q_ij - ((r_ij.^2 + x_ij.^2)).*L_ij];
Constraints = [Constraints,U_i(mpc.branch(:,1))-U_i(mpc.branch(:,2))>= -m_ij + 2*r_ij.*P_ij + 2*x_ij.*Q_ij - ((r_ij.^2 + x_ij.^2)).*L_ij];
for k=1:nlConstraints = [Constraints, cone([2*P_ij(k) 2*Q_ij(k) L_ij(k)-U_i(mpc.branch(k,1))],L_ij(k)+U_i(mpc.branch(k,1)))];
end
Constraints = [Constraints, Sij_max'.^2.*z_ij >= P_ij.^2+Q_ij.^2];%% 2.拓扑约束
Constraints = [Constraints , sum(z_ij) == nb-ns];%% 3.注入功率约束
Constraints = [Constraints, P_g>=P_g_min,P_g<=P_g_max];
Constraints = [Constraints, Q_g>=Q_g_min,Q_g<=Q_g_max];%% 4.电压上下限约束
Constraints = [Constraints, Umin <= U_i,U_i <= Umax];%% 目标函数
objective = sum(L_ij.*r_ij);%网损最小%% 设求解器
% gurobi求解器
ops=sdpsettings('verbose', 3, 'solver', 'gurobi','showprogress',1,'debug',1);
ops.gurobi.TimeLimit=600;                 % 运行时间限制为10min
ops.gurobi.MIPGap=0.01;                   % 收敛精度限制为0.01sol=optimize(Constraints,objective,ops);
value(objective)%% 分析错误标志
if sol.problem == 0disp('求解成功');
elsedisp('运行出错');yalmiperror(sol.problem)
end

数据文件如下：

function mpc = IEEE33%% MATPOWER Case Format : Version 2
mpc.version = '2';%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 10;%% bus data
%   bus_i   type    Pd  Qd  Gs  Bs  area    Vm  Va  baseKV  zone    Vmax    Vmin
mpc.bus = [ %% (Pd and Qd are specified in kW & kVAr here, converted to MW & MVAr below)1   3   0   0   0   0   1   1   0   12.66   1   1   1;2   1   100 60  0   0   1   1   0   12.66   1   1.1 0.9;3   1   90  40  0   0   1   1   0   12.66   1   1.1 0.9;4   1   120 80  0   0   1   1   0   12.66   1   1.1 0.9;5   1   60  30  0   0   1   1   0   12.66   1   1.1 0.9;6   1   60  20  0   0   1   1   0   12.66   1   1.1 0.9;7   1   200 100 0   0   1   1   0   12.66   1   1.1 0.9;8   1   200 100 0   0   1   1   0   12.66   1   1.1 0.9;9   1   60  20  0   0   1   1   0   12.66   1   1.1 0.9;10  1   60  20  0   0   1   1   0   12.66   1   1.1 0.9;11  1   45  30  0   0   1   1   0   12.66   1   1.1 0.9;12  1   60  35  0   0   1   1   0   12.66   1   1.1 0.9;13  1   60  35  0   0   1   1   0   12.66   1   1.1 0.9;14  1   120 80  0   0   1   1   0   12.66   1   1.1 0.9;15  1   60  10  0   0   1   1   0   12.66   1   1.1 0.9;16  1   60  20  0   0   1   1   0   12.66   1   1.1 0.9;17  1   60  20  0   0   1   1   0   12.66   1   1.1 0.9;18  1   90  40  0   0   1   1   0   12.66   1   1.1 0.9;19  1   90  40  0   0   1   1   0   12.66   1   1.1 0.9;20  1   90  40  0   0   1   1   0   12.66   1   1.1 0.9;21  1   90  40  0   0   1   1   0   12.66   1   1.1 0.9;22  1   90  40  0   0   1   1   0   12.66   1   1.1 0.9;23  1   90  50  0   0   1   1   0   12.66   1   1.1 0.9;24  1   420 200 0   0   1   1   0   12.66   1   1.1 0.9;25  1   420 200 0   0   1   1   0   12.66   1   1.1 0.9;26  1   60  25  0   0   1   1   0   12.66   1   1.1 0.9;27  1   60  25  0   0   1   1   0   12.66   1   1.1 0.9;28  1   60  20  0   0   1   1   0   12.66   1   1.1 0.9;29  1   120 70  0   0   1   1   0   12.66   1   1.1 0.9;30  1   200 600 0   0   1   1   0   12.66   1   1.1 0.9;31  1   150 70  0   0   1   1   0   12.66   1   1.1 0.9;32  1   210 100 0   0   1   1   0   12.66   1   1.1 0.9;33  1   60  40  0   0   1   1   0   12.66   1   1.1 0.9;
];%% generator data
%   bus Pg  Qg  Qmax    Qmin    Vg  mBase   status  Pmax    Pmin    Pc1 Pc2 Qc1min  Qc1max  Qc2min  Qc2max  ramp_agc    ramp_10 ramp_30 ramp_q  apf
mpc.gen = [1   0   0   10  -10 1   100 1   10  0   0   0   0   0   0   0   0   0   0   0   0;
];%% branch data
%   fbus    tbus    r   x   b   rateA   rateB   rateC   ratio   angle   status  angmin  angmax
mpc.branch = [  %% (r and x specified in ohms here, converted to p.u. below)1   2   0.0922  0.0470  0   0   0   0   0   0   1   -360    360;2   3   0.4930  0.2511  0   0   0   0   0   0   1   -360    360;3   4   0.3660  0.1864  0   0   0   0   0   0   1   -360    360;4   5   0.3811  0.1941  0   0   0   0   0   0   1   -360    360;5   6   0.8190  0.7070  0   0   0   0   0   0   1   -360    360;6   7   0.1872  0.6188  0   0   0   0   0   0   1   -360    360;7   8   0.7114  0.2351  0   0   0   0   0   0   1   -360    360;8   9   1.0300  0.7400  0   0   0   0   0   0   1   -360    360;9   10  1.0440  0.7400  0   0   0   0   0   0   1   -360    360;10  11  0.1966  0.0650  0   0   0   0   0   0   1   -360    360;11  12  0.3744  0.1238  0   0   0   0   0   0   1   -360    360;12  13  1.4680  1.1550  0   0   0   0   0   0   1   -360    360;13  14  0.5416  0.7129  0   0   0   0   0   0   1   -360    360;14  15  0.5910  0.5260  0   0   0   0   0   0   1   -360    360;15  16  0.7463  0.5450  0   0   0   0   0   0   1   -360    360;16  17  1.2890  1.7210  0   0   0   0   0   0   1   -360    360;17  18  0.7320  0.5740  0   0   0   0   0   0   1   -360    360;2   19  0.1640  0.1565  0   0   0   0   0   0   1   -360    360;19  20  1.5042  1.3554  0   0   0   0   0   0   1   -360    360;20  21  0.4095  0.4784  0   0   0   0   0   0   1   -360    360;21  22  0.7089  0.9373  0   0   0   0   0   0   1   -360    360;3   23  0.4512  0.3083  0   0   0   0   0   0   1   -360    360;23  24  0.8980  0.7091  0   0   0   0   0   0   1   -360    360;24  25  0.8960  0.7011  0   0   0   0   0   0   1   -360    360;6   26  0.2030  0.1034  0   0   0   0   0   0   1   -360    360;26  27  0.2842  0.1447  0   0   0   0   0   0   1   -360    360;27  28  1.0590  0.9337  0   0   0   0   0   0   1   -360    360;28  29  0.8042  0.7006  0   0   0   0   0   0   1   -360    360;29  30  0.5075  0.2585  0   0   0   0   0   0   1   -360    360;30  31  0.9744  0.9630  0   0   0   0   0   0   1   -360    360;31  32  0.3105  0.3619  0   0   0   0   0   0   1   -360    360;32  33  0.3410  0.5302  0   0   0   0   0   0   1   -360    360;21  8   2.0000  2.0000  0   0   0   0   0   0   0   -360    360;9   15  2.0000  2.0000  0   0   0   0   0   0   0   -360    360;12  22  2.0000  2.0000  0   0   0   0   0   0   0   -360    360;18  33  0.5000  0.5000  0   0   0   0   0   0   0   -360    360;25  29  0.5000  0.5000  0   0   0   0   0   0   0   -360    360;
];%%-----  OPF Data  -----%%
%% generator cost data
%   1   startup shutdown    n   x1  y1  ... xn  yn
%   2   startup shutdown    n   c(n-1)  ... c0
mpc.gencost = [2   0   0   3   0   20  0;
];%% convert branch impedances from Ohms to p.u.
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
Vbase = mpc.bus(1, BASE_KV) * 1e3;      %% in Volts
Sbase = mpc.baseMVA * 1e6;              %% in VA
mpc.branch(:, [BR_R BR_X]) = mpc.branch(:, [BR_R BR_X]) / (Vbase^2 / Sbase);%% convert loads from kW to MW
mpc.bus(:, [PD, QD]) = mpc.bus(:, [PD, QD]) / 1e3;

运行上面的代码，得到的结果如下：

提示我们优化问题无解或者是目标函数无界。按照2.4.1小节第3点所提到的，我们先把目标函数设置为空，确定一下具体是优化问题无解还是目标函数无界：

sol=optimize(Constraints,[],ops);

得到的结果如下：

OK，现在确定了是优化问题不可行而不是目标函数无界，我们可以继续调试。

1）首先只保留第一条约束，将其他约束都设为注释，查看模型是否可行：

%% 清除内存空间
clc
clear
close all
warning off
%% 系统参数
mpc = IEEE33;
nb=33;                              % 节点数
ns=1;                               % 电源节点数
nl=37;                              % 支路数
P_load=mpc.bus(:,3)/mpc.baseMVA;  % 有功负荷
Q_load=mpc.bus(:,4)/mpc.baseMVA;  % 无功负荷
r_ij=mpc.branch(:,3);              % 线路电阻
x_ij=mpc.branch(:,4);              % 线路电抗
M=1.06*1.06 - 0.94*0.94;
% 电源功率最大值
P_g_max=zeros(nb,1);
P_g_max(1)=10/mpc.baseMVA;
Q_g_max=zeros(nb,1);
Q_g_max(1)=10/mpc.baseMVA;
% 电源功率最小值
P_g_min=zeros(nb,1);
Q_g_min=[-10/mpc.baseMVA;zeros(nb-1,1)];
% 支路功率最大值
Sij_max=15/mpc.baseMVA;
% 最大电压
Umax=[1;1.1*1.1*ones(32,1)];
Umin=[1;0.95*0.95*ones(32,1)];
% 流入节点的支路
branch_to_node=zeros(nb,nl);
% 流出节点的支路
branch_from_node=zeros(nb,nl);
for k=1:nlbranch_to_node(mpc.branch(k,2),k)=1;branch_from_node(mpc.branch(k,1),k)=1;
end%% 优化变量
z_ij=binvar(nl,1);% 支路开断情况
U_i=sdpvar(nb,1);% 电压的平方
L_ij=sdpvar(nl,1);% 电流的平方
P_ij=sdpvar(nl,1);% 线路有功功率
Q_ij=sdpvar(nl,1);% 线路无功功率
P_g=sdpvar(nb,1);% 电源有功出力
Q_g=sdpvar(nb,1);% 电源无功功率%% 约束条件
Constraints = [];
%% 1.潮流约束
m_ij=(1-z_ij)*M; 
Constraints = [Constraints, P_g-P_load+branch_to_node*P_ij-branch_to_node*(L_ij.*r_ij)-branch_from_node*P_ij == 0];%% 目标函数
objective = sum(L_ij.*r_ij);%网损最小%% 设求解器
% gurobi求解器
ops=sdpsettings('verbose', 3, 'solver', 'gurobi','showprogress',1,'debug',1);
ops.gurobi.TimeLimit=600;                 % 运行时间限制为10min
ops.gurobi.MIPGap=0.01;                   % 收敛精度限制为0.01sol=optimize(Constraints, [], ops);
value(objective)%% 分析错误标志
if sol.problem == 0disp('求解成功');
elsedisp('运行出错');yalmiperror(sol.problem)
end

运行结果如下：

这个结果，说明第一条约束没问题，可以继续添加约束。

2）只保留前2条约束，将其他约束都设为注释，查看模型是否可行：

%% 清除内存空间
clc
clear
close all
warning off
%% 系统参数
mpc = IEEE33;
nb=33;                              % 节点数
ns=1;                               % 电源节点数
nl=37;                              % 支路数
P_load=mpc.bus(:,3)/mpc.baseMVA;  % 有功负荷
Q_load=mpc.bus(:,4)/mpc.baseMVA;  % 无功负荷
r_ij=mpc.branch(:,3);              % 线路电阻
x_ij=mpc.branch(:,4);              % 线路电抗
M=1.06*1.06 - 0.94*0.94;
% 电源功率最大值
P_g_max=zeros(nb,1);
P_g_max(1)=10/mpc.baseMVA;
Q_g_max=zeros(nb,1);
Q_g_max(1)=10/mpc.baseMVA;
% 电源功率最小值
P_g_min=zeros(nb,1);
Q_g_min=[-10/mpc.baseMVA;zeros(nb-1,1)];
% 支路功率最大值
Sij_max=15/mpc.baseMVA;
% 最大电压
Umax=[1;1.1*1.1*ones(32,1)];
Umin=[1;0.95*0.95*ones(32,1)];
% 流入节点的支路
branch_to_node=zeros(nb,nl);
% 流出节点的支路
branch_from_node=zeros(nb,nl);
for k=1:nlbranch_to_node(mpc.branch(k,2),k)=1;branch_from_node(mpc.branch(k,1),k)=1;
end%% 优化变量
z_ij=binvar(nl,1);% 支路开断情况
U_i=sdpvar(nb,1);% 电压的平方
L_ij=sdpvar(nl,1);% 电流的平方
P_ij=sdpvar(nl,1);% 线路有功功率
Q_ij=sdpvar(nl,1);% 线路无功功率
P_g=sdpvar(nb,1);% 电源有功出力
Q_g=sdpvar(nb,1);% 电源无功功率%% 约束条件
Constraints = [];
%% 1.潮流约束
m_ij=(1-z_ij)*M; 
Constraints = [Constraints, P_g-P_load+branch_to_node*P_ij-branch_to_node*(L_ij.*r_ij)-branch_from_node*P_ij == 0];
Constraints = [Constraints, Q_g-Q_load+branch_to_node*Q_ij-branch_to_node*(L_ij.*x_ij)-branch_from_node*Q_ij == 0];%% 目标函数
objective = sum(L_ij.*r_ij);%网损最小%% 设求解器
% gurobi求解器
ops=sdpsettings('verbose', 3, 'solver', 'gurobi','showprogress',1,'debug',1);
ops.gurobi.TimeLimit=600;                 % 运行时间限制为10min
ops.gurobi.MIPGap=0.01;                   % 收敛精度限制为0.01sol=optimize(Constraints, [], ops);
value(objective)%% 分析错误标志
if sol.problem == 0disp('求解成功');
elsedisp('运行出错');yalmiperror(sol.problem)
end

运行结果如下：

这个结果，说明前2条约束没问题且不存在矛盾，可以继续添加约束。

……

3）一直重复上述步骤，可以发现直到加入最后一条电压上下限约束之前，模型都是可行的，而加入电压上下限约束之后模型就变得不可行了。说明问题很大可能就出现在这个电压约束上。可能是电压上下限设置的太紧，我们可以尝试将节点电压标幺值的下限从0.95改成0.9，再重新运行模型。