Unbounded solutions cannot be obtained when solving LP problems that are unbounded (feasibility cut for benders decomposition)
Answered
I am now encountering a new problem that when the sub-problem is unbounded, the unbounded solution of the model cannot be called up through the solver, which leads to the failure of the model iteration (because the feasibility cut cannot be returned), I have set gurobi's InfUnbdInfo to 1 but it still does not work, I would like to consult how to get the unbounded solution of the LP problem under the unbounded.
Here are the codes:
% 加载MATPOWER数据
mpc = case39; % 假设您有一个名为 case39 的数据集
baseMVA = mpc.baseMVA; % 基准功率
busData = mpc.bus; % 节点数据
genData = mpc.gen; % 发电机数据
branchData = mpc.branch; % 支路数据
% 定义常量
F_BUS = 1;
T_BUS = 2;
BR_X = 4;
RATE_A = 6;
GEN_BUS = 1;
PMAX = 9;
PD = 3;
M = 1e8; % 大常数
z_bar=[0;0;0;0;1;0;0;0;0;0;0;0;0;1;0;0;0;0;0;1;0;0;0;0;0;0;1;1;1;0;0;1;1;1;1;1;1;1;1;0;1;0;1;1;1;1];
% 变量数量
nBranches = size(branchData, 1);
nGenerators = size(genData, 1);
nBuses = size(busData, 1);
% 创建 Gurobi 模型
model.A = sparse(0, 0);
model.obj = [];
model.rhs = [];
model.sense = '';
% 添加模型变量
mu_k_plus = 1:nBranches; % μ̅_k
mu_k_minus = nBranches + (1:nBranches); % ▁μ_k
delta_k_plus = 2 * nBranches + (1:nBranches); % δ̅_k
delta_k_minus = 3 * nBranches + (1:nBranches); % ▁δ_k
sigma_g = 4 * nBranches + (1:nGenerators); % σ̅_g
lmp_n = 4 * nBranches + nGenerators + (1:nBuses); % lmp_n
% 目标函数
model.obj = zeros(4 * nBranches + nGenerators + nBuses, 1);
% ∑_n P_(d,n) * lmp_n
model.obj(lmp_n) = busData(:, PD) / baseMVA;
% -∑_k P_k^max * (z_k)_bar * (μ̅_k + ▁μ_k)
for k = 1:nBranches
model.obj(mu_k_plus(k)) = -branchData(k, RATE_A) * z_bar(k); % P_k^max * (z_k)_bar * μ̅_k
model.obj(mu_k_minus(k)) = -branchData(k, RATE_A) * z_bar(k); % P_k^max * (z_k)_bar * ▁μ_k
end
% -∑_k (1 - (z_k)_bar) * M * (δ̅_k + ▁δ_k)
for k = 1:nBranches
model.obj(delta_k_plus(k)) = -(1 - z_bar(k)) * M; % (1 - (z_k)_bar) * M * δ̅_k
model.obj(delta_k_minus(k)) = -(1 - z_bar(k)) * M; % (1 - (z_k)_bar) * M * ▁δ_k
end
% -∑_g P_g^max * σ̅_g
for g = 1:nGenerators
model.obj(sigma_g(g)) = -genData(g, PMAX); % P_g^max * σ̅_g
end
% 添加约束
constraints = [];
rhs = [];
senses = '';
% Equation (27): 节点功率平衡
for n = 1:nBuses
delta_in = sum(delta_k_plus(branchData(:, F_BUS) == n)) - sum(delta_k_minus(branchData(:, F_BUS) == n));
delta_out = sum(delta_k_plus(branchData(:, T_BUS) == n)) - sum(delta_k_minus(branchData(:, T_BUS) == n));
constraint = zeros(1, 4 * nBranches + nGenerators + nBuses);
constraint(delta_k_plus(branchData(:, F_BUS) == n)) = 1;
constraint(delta_k_minus(branchData(:, F_BUS) == n)) = -1;
constraint(delta_k_plus(branchData(:, T_BUS) == n)) = -1;
constraint(delta_k_minus(branchData(:, T_BUS) == n)) = 1;
constraints = [constraints; constraint];
rhs = [rhs; 0];
senses = [senses; '='];
end
% Equation (28): 支路功率流
for k = 1:nBranches
fbus = branchData(k, F_BUS);
tbus = branchData(k, T_BUS);
if branchData(k, BR_X) == 0
warning('Branch %d has zero reactance, skipping.', k);
continue;
end
B_k = 1 / branchData(k, BR_X);
constraint = zeros(1, 4 * nBranches + nGenerators + nBuses);
constraint(lmp_n(fbus)) = B_k;
constraint(lmp_n(tbus)) = -B_k;
constraint(mu_k_plus(k)) = -B_k;
constraint(mu_k_minus(k)) = B_k;
constraint(delta_k_plus(k)) = -1;
constraint(delta_k_minus(k)) = 1;
constraints = [constraints; constraint];
rhs = [rhs; 0];
senses = [senses; '='];
end
% Equation (29): 发电机约束
for g = 1:nGenerators
n = genData(g, GEN_BUS);
constraint = zeros(1, 4 * nBranches + nGenerators + nBuses);
constraint(lmp_n(n)) = 1;
constraint(sigma_g(g)) = -1;
constraints = [constraints; constraint];
rhs = [rhs; genData(g, PMAX)];
senses = [senses; '<'];
end
% 非负约束
for k = 1:nBranches
% μ̅_k >= 0
constraint = zeros(1, 4 * nBranches + nGenerators + nBuses);
constraint(mu_k_plus(k)) = 1;
constraints = [constraints; constraint];
rhs = [rhs; 0];
senses = [senses; '>'];
% ▁μ_k >= 0
constraint = zeros(1, 4 * nBranches + nGenerators + nBuses);
constraint(mu_k_minus(k)) = 1;
constraints = [constraints; constraint];
rhs = [rhs; 0];
senses = [senses; '>'];
% δ̅_k >= 0
constraint = zeros(1, 4 * nBranches + nGenerators + nBuses);
constraint(delta_k_plus(k)) = 1;
constraints = [constraints; constraint];
rhs = [rhs; 0];
senses = [senses; '>'];
% ▁δ_k >= 0
constraint = zeros(1, 4 * nBranches + nGenerators + nBuses);
constraint(delta_k_minus(k)) = 1;
constraints = [constraints; constraint];
rhs = [rhs; 0];
senses = [senses; '>'];
end
for g = 1:nGenerators
% σ̅_g >= 0
constraint = zeros(1, 4 * nBranches + nGenerators + nBuses);
constraint(sigma_g(g)) = 1;
constraints = [constraints; constraint];
rhs = [rhs; 0];
senses = [senses; '>'];
end
% 将约束转化为 Gurobi 格式
model.A = sparse(constraints);
model.rhs = rhs;
model.sense = char(senses); % 转换为字符数组
% 设置求解器参数
params.outputflag = 1; % 打印输出
params.InfUnbdInfo = 1; % 查询无界射线
% 求解模型
result = gurobi(model, params);
% 检查求解状态
if strcmp(result.status, 'OPTIMAL')
fprintf('Optimal solution found.\n');
isBounded = true;
subopt = result.objval;
mu_k = result.x(mu_k_plus) + result.x(mu_k_minus); % 提取 μ̅_k 和 ▁μ_k
sigma_g = result.x(sigma_g); % 提取 σ̅_g
delta_k = result.x(delta_k_plus) + result.x(delta_k_minus); % 提取 δ̅_k 和 ▁δ_k
lmp_n = result.x(lmp_n); % 提取 lmp_n
elseif strcmp(result.status, 'INFEASIBLE')
fprintf('The problem is infeasible.\n');
subopt = inf;
isBounded = false;
mu_k = [];
sigma_g = [];
delta_k = [];
lmp_n = [];
elseif strcmp(result.status, 'UNBOUNDED')
fprintf('The problem is unbounded.\n');
isBounded = false;
% 提取无界射线
if isfield(result, 'unbddray')
unbddray = result.unbddray;
mu_k = unbddray(mu_k_plus) + unbddray(mu_k_minus);
sigma_g = unbddray(sigma_g);
delta_k = unbddray(delta_k_plus) + unbddray(delta_k_minus);
lmp_n = unbddray(lmp_n);
else
mu_k = [];
sigma_g = [];
delta_k = [];
lmp_n = [];
end
else
fprintf('Solver failed with status: %s\n', result.status);
subopt = inf;
isBounded = false;
mu_k = [];
sigma_g = [];
delta_k = [];
lmp_n = [];
end
0
-
Hi Weiqi,
Perhaps this is just the result of a "typo".
The attribute is UnbdRay, not UnbddRay.
- Riley
0
Please sign in to leave a comment.
Comments
1 comment