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，再重新运行模型。

Umin=[1;0.9*0.9*ones(32,1)];

        这时候发现可以求解出优化结果了：

        我们可以再次把目标函数添加到优化问题中，模型依旧可行。

        所以，这个优化问题不可行的原因就是电压约束设置的太紧，稍微松弛一些即可。

        对于其他不可行的优化问题，也可以按照相同的方式进行调试，找到存在问题的约束或互相矛盾的约束，再进行针对性的调整。

3.测试题

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.hqwc.cn/news/705178.html

如若内容造成侵权/违法违规/事实不符，请联系编程知识网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

中国农业大学：学硕11408复试线上涨40分，今年还会持续涨吗？中国农业大学计算机考研考情分析！

中国农业大学：学硕11408复试线上涨40分，今年还会持续涨吗？中国农业大学计算机考研考情分析！

中国农业大学&#xff08;China Agricultural University&#xff09;&#xff0c;简称“中国农大”&#xff0c;坐落于中国首都北京&#xff0c;由中华人民共和国教育部直属&#xff0c;中央直管副部级建制&#xff0c;水利部、农业部和北京市共建&#xff0c;位列国家“双一流…
阅读更多...
解决kali Linux2024无法获取动态IPv4地址（DHCP）解决方案

解决kali Linux2024无法获取动态IPv4地址（DHCP）解决方案

用root用户启动终端 进入根目录&#xff0c;选择配置文件 cd到根目录下/../etc/network找到interfaces文件 编辑interfaces文件 vi interfaces&#xff0c;编辑interfaces文件 输入如下命令 打开虚拟网络编辑器 选择虚拟机选项卡&#xff0c;编辑&#xff0c;打开虚拟网络编…
阅读更多...
Jmeter(四十一) - 从入门到精通进阶篇 - Jmeter配置文件的刨根问底 - 下篇（详解教程）

Jmeter(四十一) - 从入门到精通进阶篇 - Jmeter配置文件的刨根问底 - 下篇（详解教程）

宏哥微信粉丝群&#xff1a;https://bbs.csdn.net/topics/618423372 有兴趣的可以扫码加入 1.简介 为什么宏哥要对Jmeter的配置文件进行一下讲解了&#xff0c;因为有的童鞋或者小伙伴在测试中遇到一些需要修改配置文件的问题不是很清楚也不是很懂&#xff0c;就算修改了也是…
阅读更多...
tag-字符串：最长公共前缀

tag-字符串：最长公共前缀

题目 编写一个函数来查找字符串数组中的最长公共前缀。 如果不存在公共前缀&#xff0c;返回空字符串 “”。 示例 题解一 class Solution:def longestCommonPrefix(self, strs: List[str]) -> str:# 按照字典顺序找到strs中最大的字符串和最小的字符串str0 min(strs)st…
阅读更多...
51单片机超声波测距_液位检测_温度检测原理图PCB仿真代码

51单片机超声波测距_液位检测_温度检测原理图PCB仿真代码

目录 实物图&#xff1a; PCB ​原理图​ 仿真图 ​编辑 程序 资料下载地址&#xff1a;51单片机超声波测距-液位检测-温度检测原理图PCB仿真代码 主控为stc89c52,通过ds18b20进行温度采集&#xff0c;超声波测距&#xff0c;距离不可以超过1m&#xff0c;通过按键可以设…
阅读更多...
Java面试八股之List、Set、Map和Queue之间的区别

Java面试八股之List、Set、Map和Queue之间的区别

Java中List、Set、Map和Queue之间的区别 List List 是一个有序且允许重复元素的集合。它提供了按索引访问元素的能力&#xff0c;也就是说&#xff0c;你可以通过元素的插入位置或指定的索引来精确地访问、添加或删除元素。List 的典型实现包括 ArrayList 和 LinkedList。 特…
阅读更多...
喜茶·茶坊黑金首店入驻北京三里屯，率先引入珍稀娟姗奶制茶

喜茶·茶坊黑金首店入驻北京三里屯，率先引入珍稀娟姗奶制茶

发布 | 大力财经 近日&#xff0c;喜茶茶坊 BLACK 在北京三里屯开业&#xff0c;这是喜茶新业态的首家黑金店型。该店在延续喜茶茶坊“鲜、茶、纯”的精品茗茶特色和宋代茶文化审美意趣的基础上&#xff0c;首次升级呈现了铜锅手煮烹茶工艺、娟姗牛乳制茶等创新尝试&#xff0…
阅读更多...
【C++ 高阶数据结构 Test】AVL ~ 二叉搜索树

【C++ 高阶数据结构 Test】AVL ~ 二叉搜索树

文章目录 1. AVL 树概念2. AVL 树节点的定义3. AVL树的插入4. AVL树的旋转4.1 新节点插入较高左子树的左侧---左左&#xff1a;右单旋4.2 新节点插入较高右子树的右侧---右右&#xff1a;左单旋4.3 新节点插入较高左子树的右侧---左右&#xff1a;先左单旋再右单旋4.4 新节点插…
阅读更多...
电工能混到这份上

电工能混到这份上

最近看到某电工师傅发了一篇帖子&#xff0c;大致内容是他在处理一个简单故障的时候居然花了很长的时间。我们一起来看看他遇到的是什么故障吧! plc 控制的一台设备&#xff0c;行走部分靠 2 个脚踏开关控制&#xff08;内部开关量控制方向&#xff0c;电位器控制速度&#xff…
阅读更多...
微信小程序快速开发-基础内容（内容真的又多又干货）

微信小程序快速开发-基础内容（内容真的又多又干货）

目录 实现横向布局效果 实现滚动效果 实现轮播图效果 实现文本长按选中复制效果 渲染 HTML 标签 按钮组件的使用效果 图片组件的使用效果 Mustache 语法 动态绑定内容&#xff08;定义变量&#xff0c;渲染变量&#xff09; 动态绑定属性&#xff08;将属性定义为变量…
阅读更多...
用SwitchHosts模拟本地域名解析访问

用SwitchHosts模拟本地域名解析访问

一.用SwitchHosts模拟本地域名解析访问 1.下载地址 https://download.csdn.net/download/jinhuding/89313168 2.使用截图
阅读更多...
软考-软件工程

软考-软件工程

软件工程概述 软件工程指的是应用计算机科学、数学及管理科学等原理&#xff0c;以工程化的原则和方法来解决软件 问题的工程&#xff0c;目的是提高软件生产率、提高软件质量、降低软件成本。 概述&#xff1a; 软件开发模型&#xff1a;指导软件开发的体系 需求分析确定软件…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023