Numerical trouble encountered when solving convex QCQP
回答済み
Dear Sir/Madam,
I encountered a numerical problem while solving a convex QCQP problem with gurobi in Python.
I am working on a difference-of-convex algorithm. When the hyperparameters of iteration penalty coefficient \rho is set to a bit larger value and the truncation hyperparameter \delta is set to a bit smaller value, then the error will occur after several iterations:
"Barrier performed 52 iterations in 0.02 seconds (0.00 work units)
Numerical trouble encountered
Model may be infeasible or unbounded. Consider using the
homogeneous algorithm (through parameter 'BarHomogeneous')"
Core code is here. When \delta and \rho are set to values like a=0.05, b=10, then the model can be solved appropriately. However, when a=0.02, b=15, then the above NUMERICAL status code is encountered.
Looking forward to your advice. Thank you.
class new_DCA:
def __init__(self, delta, rho, llambda, T, m, p, q, Delta, K, c, Cbar):
self.delta = delta
self.rho = rho
self.llambda = llambda
self.T = T
self.m = m
self.p = p
self.q = q
self.Delta = Delta
self.I = len(m)
self.K = K
self.c = c
self.Cbar = Cbar
def truncation_index(self, x, delta_now):
if delta_now - x >= 0:
return 1
else:
return 0
def main_algorithm(self, X):
optimal = []
for k in range(len(self.delta)):
delta_now = self.delta[k]
model = gp.Model()
X_var = model.addVars(self.I, self.T+1, lb=0,ub=1, vtype=GRB.CONTINUOUS)
penalty = model.addVar()
model.addLConstr(penalty >= gp.quicksum([self.c[i] * gp.quicksum([1-(1-X_var[i,t]/delta_now)*self.truncation_index(X[i,t], delta_now=delta_now) for t in range(1,self.T+1)])
+ self.K[i]*(1-(1-X_var[i,self.T]/delta_now)*self.truncation_index(X[i,self.T], delta_now=delta_now)) for i in range(self.I)]) - self.Cbar)
for i in range(self.I):
model.addLConstr(X_var[i,0]==0)
for t in range(1, self.T+1):
model.addLConstr(X_var[i, t-1] <= X_var[i, t])
model.addConstr(X_var[i,t] <= X_var[i, t-1]*(1 - p[i])
+ (p[i] + q[i]*gp.quicksum([self.Delta[i,j]*X_var[j, t-1] for j in range(self.I)]))
- q[i]*self.Delta[i,i]*(0.5*gp.quicksum([self.Delta[i,j]/self.Delta[i,i]*X_var[j, t-1] if j != i else 0 for j in range(self.I)]) + X_var[i, t-1]) * (0.5*gp.quicksum([self.Delta[i,j]/self.Delta[i,i]*X_var[j, t-1] if j != i else 0 for j in range(self.I)]) + X_var[i, t-1])
- 0.25*q[i]/self.Delta[i,i]*sum([self.Delta[i,j]*X[j,t-1] if j != i else 0 for j in range(self.I)]) * sum([self.Delta[i,j]*X[j,t-1] if j != i else 0 for j in range(self.I)])
+ 0.5*q[i]/self.Delta[i,i]*sum([self.Delta[i,j]*X[j,t-1] if j != i else 0 for j in range(self.I)]) * gp.quicksum([self.Delta[i,j]*X_var[j, t-1] if j != i else 0 for j in range(self.I)]))
model.setObjective(-1 * gp.quicksum([self.m[i] * X_var[i, T] for i in range(self.I)]) + 1/2/self.llambda * gp.quicksum([(X_var[i,t] - X[i,t])**2 for i in range(self.I) for t in range(1,self.T+1)]) + self.rho[k] * penalty)
model.update()
model.optimize()
if model.status == GRB.OPTIMAL:
for i in range(self.I):
for t in range(self.T+1):
X[i,t]=X_var[i,t].x
obj_value = -1 * sum([self.m[i] * X[i, T] for i in range(self.I)])
optimal.append(round(obj_value,4))
else:
return False
return optimal, X
if __name__ == '__main__':
m = [1000 for i in range(4)]
p = [0.02 for i in range(4)]
q = [0.4 for i in range(4)]
Delta = np.array([[0.25,0.25,0.25,0.25],
[0.25,0.25,0.25,0.25],
[0.25,0.25,0.25,0.25],
[0.25,0.25,0.25,0.25]])
T = 3
c = [1,1,1,1]
K = [2,1,1,1]
Cbar = 6.5
iteration = 30
a = 0.02
b = 15
deltas = [a*0.9**i for i in range(iteration)]
rhos = [b*i**1.2 for i in range(1,iteration+1)]
llambda = 100
X = np.array([[0.0 for t in range(T+1)] for i in range(4)])
examp = new_DCA(deltas, rhos, llambda, T, m, p ,q, Delta, K, c, Cbar)
result, Xs = examp.main_algorithm(deepcopy(X))
-
正式なコメント
-
Gurobi Staff
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum, or try Gurobot, our chatbot interface offering instant, expert-level support. -
-
-
Hi Braun,
Your objective function range increases with every iteration until it gets too large to be manageable. At the last iteration, we get:Coefficient statistics:
Matrix range [1e+00, 3e+03]
QMatrix range [3e-02, 1e-01]
QLMatrix range [1e-01, 1e+00]
Objective range [2e-14, 1e+03]
QObjective range [1e-02, 1e-02]
Bounds range [1e+00, 1e+00]
RHS range [5e-01, 5e-01]
QRHS range [2e-02, 2e-02]Changing lambda to 10 also works.
You need to try and scale these coefficients to reduce this range.
Please see: Guidelines for Numerical Issues
As suggested in that guide, using Presolve = 0 also works for this model.
Cheers,
David0 -
投稿コメントは受け付けていません。
コメント
2件のコメント