Skip to main content

Numerical trouble encountered when solving convex QCQP

Answered
Braun Reiner
  • Gurobi-versary
  • First Question

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))
0

Comments

2 comments

Date Votes
  • Official comment
    Simranjit Kaur
    • 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.
  • David Torres Sanchez
    • Gurobi Staff
    Edited

    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, 
    David

    0

Post is closed for comments.

Didn't find what you were looking for?

New post