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.xdnf.cn/news/1424467.html

如若内容造成侵权/违法违规/事实不符，请联系一条长河网进行投诉反馈，一经查实，立即删除！

相关文章

Tomcat无法连通的调试方法1-service方式无法连通

Tomcat无法连通的调试方法1-service方式无法连通

作者&#xff1a;私语茶馆 1.局域网Tomcat服务不通 组网如下&#xff1a; 问题&#xff1a; Tomcat Server 服务方式启动后&#xff0c;无法访问&#xff0c;但命令行方式启动可以。IP地址都在同网段或不同网段现象都一样。 2.Tomcat 服务安装与调试 在Windows下&#xff0c;…
阅读更多...
Cadence 16.6 绘制PCB封装时总是卡死的解决方法

Cadence 16.6 绘制PCB封装时总是卡死的解决方法

Cadence 16.6 绘制PCB封装时总是卡死的解决方法 在用Cadence 16.6 PCB Editor绘制PCB封装时候&#xff0c;绘制一步卡死一步&#xff0c;不知道怎么回事儿&#xff0c;在咨询公司IT后&#xff0c;发现是WIN系统自带输入法的某些热键与PCB Editor有冲突&#xff0c;导致卡死。 …
阅读更多...
第十四届蓝桥杯大赛软件赛国赛C/C++ 大学 B 组 拼数字

第十四届蓝桥杯大赛软件赛国赛C/C++ 大学 B 组 拼数字

//bfs只能过40%。 #include<bits/stdc.h> using namespace std; #define int long long int a,b,c,dp[2028]; struct s {int x,y,z;string m; }; map<vector<int>,int>k; signed main() {ios::sync_with_stdio(false);cin.tie(0),cout.tie(0);cin>>a…
阅读更多...
【SAP ME 36】NetWearver数据源重定向（Oracle）

【SAP ME 36】NetWearver数据源重定向（Oracle）

禁止废话，直接上图！！！
阅读更多...
[链表专题]力扣141, 142

[链表专题]力扣141, 142

1. 力扣141 : 环形链表 题 : 给你一个链表的头节点 head &#xff0c;判断链表中是否有环。 如果链表中有某个节点&#xff0c;可以通过连续跟踪 next 指针再次到达&#xff0c;则链表中存在环。 为了表示给定链表中的环&#xff0c;评测系统内部使用整数 pos 来表示链表尾…
阅读更多...
git 基础

git 基础

【一】介绍 &#xff08;1&#xff09;软件开发模式 瀑布式开发 线性、顺序的开发流程&#xff0c;严格按照需求分析、设计、编码、测试、交付等阶段进行。每个阶段都必须在前一个阶段完成后才能开始&#xff0c;因此开发周期长&#xff0c;难以应对需求变更。适用于需求明确…
阅读更多...
2024网上可申请离婚，无需对方同意！

2024网上可申请离婚，无需对方同意！

&#x1f383;很多客户决定离婚之后却因为不了解离婚流程没准备好所需材料&#xff0c;导致离婚失败&#xff0c;或者无故被对方e意拖延&#xff0c;无计可施&#xff0c;无可奈何&#xff01; &#x1f383;别怕&#xff0c;2024年离婚新规定已发布&#xff0c;离婚变的简单了…
阅读更多...
利用管道通信(pipe)测量进程间的上下文切换(context switch)开销

利用管道通信(pipe)测量进程间的上下文切换(context switch)开销

利用管道通信(pipe)测量进程间的上下文切换(context switch)开销 《https://pages.cs.wisc.edu/~remzi/OSTEP/cpu-mechanisms.pdf》 Measuring the cost of a context switch is a little trickier. The lmbench benchmark does so by running two processes on a single CPU…
阅读更多...
山东大学计算机考研数据分析，初复试占比6：4，复试内容不少得花精力准备！

山东大学计算机考研数据分析，初复试占比6：4，复试内容不少得花精力准备！

山东大学&#xff08;ShandongUniversity&#xff09;&#xff0c;简称山大&#xff0c;位于中国山东&#xff0c;是中华人民共和国教育部直属的综合性全国重点大学&#xff0c;是国家“211工程”、“985工程”重点建设院校&#xff0c;入选“111计划”、“珠峰计划”、“卓越工…
阅读更多...
kafka用java收发消息

kafka用java收发消息

用java客户端代码来对kafka收发消息 具体代码如下 package com.cool.interesting.kafka;import org.apache.kafka.clients.consumer.ConsumerConfig; import org.apache.kafka.clients.consumer.ConsumerRecord; import org.apache.kafka.clients.consumer.ConsumerRecords; i…
阅读更多...
基于springboot+vue+Mysql的大学生社团活动平台

基于springboot+vue+Mysql的大学生社团活动平台

开发语言&#xff1a;Java框架&#xff1a;springbootJDK版本&#xff1a;JDK1.8服务器&#xff1a;tomcat7数据库&#xff1a;mysql 5.7&#xff08;一定要5.7版本&#xff09;数据库工具&#xff1a;Navicat11开发软件&#xff1a;eclipse/myeclipse/ideaMaven包&#xff1a;…
阅读更多...
电脑缺失api-ms-win-crt-runtime-l1-1-0.dll文件的几种修复方法

电脑缺失api-ms-win-crt-runtime-l1-1-0.dll文件的几种修复方法

当您在使用电脑过程中遇到程序启动失败&#xff0c;提示缺少“api-ms-win-crt-runtime-l1-1-0.dll”文件时&#xff0c;不必过于焦虑&#xff0c;此问题通常与Windows系统的Visual C Redistributable组件未正确安装或损坏有关。小编将介绍5种修复电脑缺失api-ms-win-crt-runtim…
阅读更多...
STM32-09-IWDG

STM32-09-IWDG

文章目录 STM32 IWDG1. IWDG2. IWDG框图3. IWDG寄存器4. IWDG寄存器操作步骤5. IWDG溢出时间计算6. IWDG配置步骤7. 代码实现 STM32 IWDG 1. IWDG IWDG Independent watchdog&#xff0c;即独立看门狗&#xff0c;本质上是一个定时器&#xff0c;这个定时器有一个输出端&#…
阅读更多...
webpack生成模块关系依赖图示例：查看构建产物的组成部分 依赖关系图

webpack生成模块关系依赖图示例：查看构建产物的组成部分 依赖关系图

npm i -D webpack-bundle-analyzer core-js babel-loaderwebpack.config.js const BundleAnalyzerPlugin require(webpack-bundle-analyzer).BundleAnalyzerPlugin; module.exports {entry: ./src/index.js,output: {filename: main.js,},// mode: production, // 或者 produ…
阅读更多...
学前端网络安全这块还不懂？细说CSRF

学前端网络安全这块还不懂？细说CSRF

什么是CSRF&#xff1f; 举个栗子&#xff0c;比如我们需要在某个博客上删除一个文章&#xff0c;攻击者首先在自己的域构造一个页面&#xff0c;使用了一个img标签&#xff0c;其地址指向了删除博客的链接。攻击者诱使目标用户&#xff0c;也就是博客主访问这个页面&#xff…
阅读更多...
如何利用命令提示符列出文件？这里提供了几个实例供你参考

如何利用命令提示符列出文件？这里提供了几个实例供你参考

序言 什么命令可以用来列出目录中的文件&#xff1f;如何在命令提示符Windows 10/11中列出文件&#xff1f;很多人对这些问题感到困惑。在这篇文章中&#xff0c;我们详细解释了命令提示符列出文件的主题。 CMD&#xff08;命令提示符&#xff09;是一个功能强大的Windows内置…
阅读更多...
Android实践：查看Activity信息

Android实践：查看Activity信息

问题&#xff1a;本地Android SDK的monitor无法正常运行&#xff0c;看不了进程相关信息&#xff0c;确认当前显示Activity十分不便 解决办法&#xff1a;使用adb shell指令可以快速查看 命令&#xff1a; adb shell dumpsys activity activities 这个命令用于获取Android设…
阅读更多...
8个迹象表明你需要一台新笔记本电脑，看一下你的笔记本是否有其中一个

8个迹象表明你需要一台新笔记本电脑，看一下你的笔记本是否有其中一个

序言 当你第一次打开你的笔记本电脑的盒子时，它会以最高性能运行，电池寿命更长，过热最小，资源使用效率高。然而，随着笔记本电脑的老化，它将不能满足预期用途。以下几个迹象表明，可能是时候寻找并投资一款新设备了。 你的设备不再具有预期用途 如果你的笔记本电脑不再…
阅读更多...
大模型日报2024-05-15

大模型日报2024-05-15

大模型日报 2024-05-15 大模型资讯 OpenAI推出全新AI模型GPT-4o&#xff0c;具备文本、图像和音频处理能力 摘要: OpenAI公司继ChatGPT后&#xff0c;最新推出了名为GPT-4o的AI模型。这一模型不仅能够理解和生成文本&#xff0c;还新增了图像和音频的解释及生成功能。GPT-4o作为…
阅读更多...
战网国际服怎么下载 暴雪战网一键下载安装图文教程

战网国际服怎么下载 暴雪战网一键下载安装图文教程

战网国际版&#xff0c;或称为Battle.net全球版&#xff0c;是暴雪娱乐构建的一项跨越国界的综合游戏交流平台&#xff0c;它无视地理限制&#xff0c;旨在服务全球每一个角落的游戏爱好者。不同于地区专属版本&#xff0c;国际版为玩家开启了一扇无门槛的大门&#xff0c;让每…
阅读更多...
最新文章

Copyright @ 2022~2023