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.mzph.cn/web/12910.shtml

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

相关文章

PS滤镜插件Camera Raw 15.4升级，开启智能修图

PS滤镜插件Camera Raw 15.4升级，开启智能修图

前段时间Adobe 更新了photoshop 的智能AI填充功能&#xff0c;深受很多设计师朋友的喜爱。Camera Raw作为PS的一个滤镜插件对RAW图片的处理上面有一定的优势&#xff0c;Camera Raw 15.4升级了&#xff0c;开启智能修图木事&#xff0c;一起来看看吧&#xff01; Camera Raw滤镜…
阅读更多...
【2024华为HCIP831 | 高级网络工程师之路】刷题日记（18）

【2024华为HCIP831 | 高级网络工程师之路】刷题日记（18）

个人名片&#xff1a;&#x1faaa; &#x1f43c;作者简介&#xff1a;一名大三在校生&#xff0c;喜欢AI编程&#x1f38b; &#x1f43b;‍❄️个人主页&#x1f947;&#xff1a;落798. &#x1f43c;个人WeChat&#xff1a;hmmwx53 &#x1f54a;️系列专栏&#xff1a;&a…
阅读更多...
ClassificationPrimitive 内部原理

ClassificationPrimitive 内部原理

ClassificationPrimitive 内部原理 发明 ClassificationPrimitive的真是个天才。其原理是利用 webgl 的模板缓冲区实现。 渲染两次, 首先是绘制模板, 然后绘制真正的内容。 示意图: function createClass() {const { program, uniforms } WebGLProgram.buildPrograms(gl, …
阅读更多...
代码随想录算法训练营第36期DAY22

代码随想录算法训练营第36期DAY22

DAY22 654最大二叉树 自己做的时候忽略了&#xff1a;nums.length>1的题给条件。所以每次递归都要判断是否size()>1&#xff0c;不要空的。 /** * Definition for a binary tree node. * struct TreeNode { * int val; * TreeNode *left; * TreeNode *rig…
阅读更多...
牛客网刷题 | BC84 牛牛学数列2

牛客网刷题 | BC84 牛牛学数列2

目前主要分为三个专栏&#xff0c;后续还会添加&#xff1a; 专栏如下&#xff1a; C语言刷题解析 C语言系列文章 我的成长经历 感谢阅读&#xff01; 初来乍到&#xff0c;如有错误请指出&#xff0c;感谢&#xff01; 描述 这次牛牛又换了个数…
阅读更多...
「AIGC」Python实现tokens算法

「AIGC」Python实现tokens算法

本文主要介绍通过python实现tokens统计，避免重复调用openai等官方api，开源节流。 一、设计思路 初始化tokenizer使用tokenizer将文本转换为tokens计算token的数量二、业务场景 2.1 首次加载依赖 2.2 执行业务逻辑 三、核心代码 from transformers import AutoTokenizer imp…
阅读更多...
React: memo

React: memo

React.memo 允许你的组件在 props 没有改变的情况下跳过重新渲染。 const MemoizedComponent memo(SomeComponent, arePropsEqual?)React 通常在其父组件重新渲染时重新渲染一个组件。你可以使用 memo 创建一个组件&#xff0c;当它的父组件重新渲染时&#xff0c;只要它的新…
阅读更多...
《海峡科技与产业》是什么级别的期刊？是正规期刊吗？能评职称吗？

《海峡科技与产业》是什么级别的期刊？是正规期刊吗？能评职称吗？

问题解答 问&#xff1a;《海峡科技与产业》期刊是什么级别&#xff1f; 答&#xff1a;国家级 主管单位&#xff1a;中华人民共和国科学技术部 主办单位&#xff1a;科技部海峡两岸科学技术交流中心 问&#xff1a;《海峡科技与产业》影响因子&#xff1f; 答&#xff1a;…
阅读更多...
随笔：棋友们

随笔：棋友们

我是在小学二年级学会中国象棋的&#xff0c;准确说&#xff0c;是学会象棋的下棋规则的&#xff0c;师傅是二舅。我最早的对手就是同学波仔。波仔比我略早学会象棋&#xff0c;总用连珠炮欺负我&#xff0c;开局几步棋就把我将死。我不知道怎么破解。轮到我先走时&#xff0c;…
阅读更多...
扭亏为盈的赛力斯，真正进入稳态了吗？

扭亏为盈的赛力斯，真正进入稳态了吗？

“72小时内大定破1万台”。5月15日&#xff0c;问界新M5开启全国大规模交付&#xff0c;从当前取得的成绩来看&#xff0c;赛力斯的“富贵”似乎还将延续。 其实&#xff0c;此前基于问界新M7等车型的爆火&#xff0c;赛力斯已经找到了创收轨道。财报显示&#xff0c;2024年一…
阅读更多...
alist网盘自动同步

alist网盘自动同步

alist网盘自动同步 alist可以设置目录定时转存到各个网盘&#xff0c;做到夸网盘&#xff0c;多备份的效果可以将自己挂载的alist 下的各个目录相互间进行同步&#xff0c;原理是采用alist原始api调用执行同步原理1.匹配文件名称是否相同,2.文件大小是否相同&#xff0c;相同会…
阅读更多...
一文详细解析Google编码规范工具cpplint的下载安装与使用

一文详细解析Google编码规范工具cpplint的下载安装与使用

目录 一、什么是cpplint 二、cpplint能实现的功能 三、cpplint的下载与使用 1、配置python环境 2、安装cpplint 四、cpplint常用命令讲解 1、常用命令查看 2、常用命令详解 3、命令使用方式 五、 cpplint的实用技巧 1、集成cpplint 1.1、修改调用接口. 1.2、直接把…
阅读更多...
数据结构（C）：树的概念和二叉树初见

数据结构（C）：树的概念和二叉树初见

目录 &#x1f37a;0.前言 1.树概念及结构 2.认识一棵树 3.树的表示 3.1树在实际中的运用&#xff08;表示文件系统的目录树结构&#xff09; 4.二叉树 4.1特殊的二叉树 4.2二叉树的性质 &#x1f48e;5.结束语 &#x1f37a;0.前言 言C之言&#xff0c;聊C之识&…
阅读更多...
卷积模型的剪枝、蒸馏---蒸馏篇--NST特征蒸馏（以deeplabv3+为例）

卷积模型的剪枝、蒸馏---蒸馏篇--NST特征蒸馏（以deeplabv3+为例）

本文使用NST特征蒸馏实现deeplabv3+模型对剪枝后模型的蒸馏过程； 一、NST特征蒸馏简介 下面是两张叠加了热力图（heat map）的图片，从图中很容易看出这两个神经元具有很强的选择性：左图的神经元对猴子的脸部非常敏感，右侧的神经元对字符非常敏感。这种激活实际上意味着神经…
阅读更多...
C++ 日志库 log4cpp 编译、压测及其范例代码 [全流程手工实践]

C++ 日志库 log4cpp 编译、压测及其范例代码 [全流程手工实践]

文章目录 一、 log4cpp官网二、下载三、编译1.目录结构如下2.configure 编译3.cmake 编译 四、测试五、压测源码及结果1.运行环境信息2.压测源码3.压测结果 文章内容&#xff1a;包含了对其linux上的完整使用流程&#xff0c;下载、编译、安装、测试用例尝试、以及一份自己写好…
阅读更多...
IT革新狂潮：引领未来的技术趋势

IT革新狂潮：引领未来的技术趋势

方向一&#xff1a;技术革新与行业应用 当前现状&#xff1a; 量子计算&#xff1a;量子计算的研究正在加速&#xff0c;尽管目前仍处于初级阶段&#xff0c;但其在药物研发、加密技术和材料科学等领域的应用潜力已被广泛认可。 虚拟现实&#xff08;VR&#xff09;与增强现实…
阅读更多...
湖南大学OS-2018期末考试（不含解析）

湖南大学OS-2018期末考试（不含解析）

前言 不知道哪里翻出来的一张&#xff0c;看着确实像期末考卷&#xff0c;暂且放一下。或许做过&#xff0c;或许没做过。 总之答案不记得了。做完可以评论区发一下或者找我发出来。 共6道大题。 一、(30%) 1. &#xff08;6%&#xff09; 进程间通信的两种方法分别是什么&…
阅读更多...
DRF 跨域问题

DRF 跨域问题

【一】说明 CORS&#xff08;跨来源资源共享&#xff0c;Cross-Origin Resource Sharing&#xff09;是一种浏览器技术的规范&#xff0c;旨在解决浏览器同源策略&#xff08;Same-Origin Policy&#xff09;的限制&#xff0c;使得Web服务可以从不同的网域&#xff08;源&…
阅读更多...
error Error: certificate has expired

error Error: certificate has expired

用yarn命令安装依赖的时候遇到报错&#xff1a; 原因&#xff1a;可能是开了服务器代理访问导致ssl安全证书失效 解决方法&#xff1a; 在终端输入 yarn config set "strict-ssl" false -g yarn config set "strict-ssl" false -g 然后再安装依赖就不…
阅读更多...
RS2227XN功能和参数介绍及PDF资料

RS2227XN功能和参数介绍及PDF资料

RS2227XN是一款模拟开关/多路复用器 品牌: RUNIC(润石) 封装: MSOP-10 描述: USB2.0高速模拟开关 开关电路: 双刀双掷(DPDT) 通道数: 2 工作电压: 1.8V~5.5V 导通电阻(RonVCC): 10Ω 功能&#xff1a;模拟开关/多路复用器 USB2.0高速模拟开关 工作电压范围&#xff1a;1.8V ~ 5…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023