Numerical trouble encountered when solving convex QCQP

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