Cvx_status=unbounded in Benders subproblem

The Benders Subproblem is as follows:

1590684673(1).png830×251 8.69 KB

Soving the problem by Python+gurobi ，the output is as follows：

Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (win64)
Optimize a model with 45 rows, 58 columns and 126 nonzeros
Model fingerprint: 0xc76df5b0
Coefficient statistics:
Matrix range [9e-01, 2e+02]
Objective range [1e+00, 2e+03]
Bounds range [0e+00, 0e+00]
RHS range [5e-01, 8e-01]
Presolve removed 45 rows and 58 columns
Presolve time: 0.01s
Presolve: All rows and columns removed
Iteration Objective Primal Inf. Dual Inf. Time
0 1.2182075e+04 0.000000e+00 0.000000e+00 0s

Solved in 0 iterations and 0.01 seconds
Optimal objective 1.218207529e+04

Where Model.satus=2,which means the optimal solution is available. The optimal objective is equal to 12182.The problem is bounded.

Soving the problem by Matlab+cvx+mosek ，the output is as follows：
Calling Mosek 9.1.9: 115 variables, 50 equality constraints

MOSEK Version 9.1.9 (Build date: 2019-11-21 11:34:40)
Copyright © MOSEK ApS, Denmark. WWW: mosek.com
Platform: Windows/64-X86

Problem
Name :
Objective sense : min
Type : LO (linear optimization problem)
Constraints : 50
Cones : 0
Scalar variables : 115
Matrix variables : 0
Integer variables : 0

Optimizer started.
Presolve started.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 1 time : 0.00
Lin. dep. - tries : 0 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.03
Optimizer terminated. Time: 0.13

Interior-point solution summary
Problem status : DUAL_INFEASIBLE
Solution status : DUAL_INFEASIBLE_CER
Primal. obj: -1.6005008384e+03 nrm: 2e+02 Viol. con: 0e+00 var: 0e+00

Basic solution summary
Problem status : DUAL_INFEASIBLE
Solution status : DUAL_INFEASIBLE_CER
Primal. obj: -1.6005008384e+03 nrm: 2e+02 Viol. con: 0e+00 var: 0e+00
Optimizer summary
Optimizer - time: 0.13
Interior-point - iterations : 0 time: 0.05
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00

Status: Unbounded
Optimal value (cvx_optval): +Inf

Where the cvx_status=unbounded, which means the optimal solution doesn’t exist.

Question: The results obtained by the two solvers are different. And we guess whether it is caused by cvx?

We don’t even know that your two codes are mathematically equivalent, leaving aside any transformations which might be done in the different environments.

Thank you for answering. The two codes are as follows.
#python+gurobi
from numpy import*
import numpy as np
from gurobipy import*
import pandas as pd
#create a new model
model = Model(“Benders_SP”)
#Initialize coefficient matrix
model.Params.InfUnbdInfo = 1
excel_path= r’D:\software\gurobi_9\win64\examples\python\Benders\data.xlsx’
x=pd.read_excel(excel_path,‘x’)
y=pd.read_excel(excel_path,‘y’)
A=pd.read_excel(excel_path,‘A’)
Cc=0.1
Cg=[0.746547776,0.755680707,0.603636527,0.491336073,0.770983328,0.480292946,0.668648765,0.588128497,0.737004972,0.41018244]
d=pd.read_excel(excel_path,‘d’)
d_row=[1921.200755,1839.120405,1989.963025,1732.989174,1803.393696,1997.270359,1933.796329,1853.403336,1855.891279,1950.227868]
GAMA=np.array([[0,0,0,0,0]]).T
p=0.09
Q=[550,550,550,550,550,550,550,550,550,550]
rc=7.5
rd=6.25
xitao=0.5
gammaC=0.95
gammaD=0.9
gammaZ=0.9
#create variables
alpha={}
beta={}
m={}
r={}
for t in range(0,9):
alpha[t]=model.addVar(name=“alpha(%s)”%(t))
beta[t]=model.addVar(lb=0,name=“beta(%s)”%(t))
for i in range(0,4):
m[i]=model.addVar(lb=0,name=“m(%s)”%(i))
for i in range(0,4):
for t in range(0,9):
r[i,t]=model.addVar(lb=0,name=“r(%s,%s)”%(i,t))

#create constraints
model.addConstrs(gammaZ*alpha[t]-beta[t] == Cg[t] for t in range(0,9))
model.addConstrs(m[i]+r[i,t] >= alpha[t]*d.iloc[i,t]*xitao for i in range(0,4) for t in range(0,9))

#set objective
obj = LinExpr()
obj += sum(alpha[t](gammaCsum(x.iloc[i,t]rc for i in range(0,4))-gammaDsum(y.iloc[i,t]*rd for i in range(0,4))) for t in range(0,9))
obj += sum(alpha[t]*d_row[t] for t in range(0,9))
obj += sum(-m[i]*GAMA[i] for i in range(0,4))
obj += sum(-beta[t]*Q[t] for t in range(0,9))
obj += sum(r[i,t] for i in range(0,4) for t in range(0,9))
model.setObjective(obj, GRB.MINIMIZE)
model.optimize()
a=model.getAttr(GRB.Attr.UnbdRay, model.getVars())
print(a)

#matlab+cvx+mosek

function [alpha1,beta1,r,m,phi,BD_status,UB1]=Benders_SP(x,y, A, Aa, B, Bt, Cc, Cg, cmu, cxigma, d, GAMA, i, j, p, pxita, Q, rc, rd, t, t3, v, v3, w, xitao, gammaC, gammaD, gammaZ)
parameters;
cvx_begin
cvx_solver mosek
%create variables————————————————
variable alpha1(1,t)
variable beta1(1,t) nonnegative
variable phi(1,1) nonnegative
variable m(i,1) nonnegative
variable r(i,t) nonnegative
%set objective function——————————————————
maximize (sum(alpha1.(gammaCsum(x.rc,1)-gammaDsum(y.rd,1)))+sum(alpha1.sum(d,1))-sum(m.GAMA+sum(r,2))-sum(beta1.Q))
% (sum(alpha1.(gammaCsum(rcx,1)-gammaDsum(rdy,1)))-sum(beta1.Q)+phi)
subject to
%Constraint(4.24)、(4.27)——————————————————
for t1=1:t
gammaZalpha1(t1)-beta1(t1) == Cg(t1);
end
% gammaZ.(alpha1)-beta1 == Cg;
% sum(alpha1.*sum(d,1))-sum(m.*GAMA+sum(r,2))>=phi;
for t1=1:t
for i1=1:i
m(i1)+r(i1,t1)>=alpha1(t1)*d(i1,t1).*xitao;
end
end

cvx_end
BD_status=cvx_status;
UB1=cvx_optval;function [alpha1,beta1,r,m,phi,BD_status,UB1]=Benders_SP(x,y, A, Aa, B, Bt, Cc, Cg, cmu, cxigma, d, GAMA, i, j, p, pxita, Q, rc, rd, t, t3, v, v3, w, xitao, gammaC, gammaD, gammaZ)
parameters;
cvx_begin
cvx_solver mosek
%create variables————————————————
variable alpha1(1,t)
variable beta1(1,t) nonnegative
variable phi(1,1) nonnegative
variable m(i,1) nonnegative
variable r(i,t) nonnegative
%set objective function——————————————————
maximize (sum(alpha1.(gammaCsum(x.rc,1)-gammaDsum(y.rd,1)))+sum(alpha1.sum(d,1))-sum(m.GAMA+sum(r,2))-sum(beta1.Q))
% (sum(alpha1.(gammaCsum(rcx,1)-gammaDsum(rdy,1)))-sum(beta1.Q)+phi)
subject to
%Constraint(4.24)、(4.27)——————————————————
for t1=1:t
gammaZalpha1(t1)-beta1(t1) == Cg(t1);
end
% gammaZ.(alpha1)-beta1 == Cg;
% sum(alpha1.*sum(d,1))-sum(m.*GAMA+sum(r,2))>=phi;
for t1=1:t
for i1=1:i
m(i1)+r(i1,t1)>=alpha1(t1)*d(i1,t1).*xitao;
end
end

cvx_end
BD_status=cvx_status;
UB1=cvx_optval;

Based on the output I am 99.9999% sure that the problem Mosek is fed is dual infeasible.