Skip to main content

Gurobi cannot find the solution after one day

Answered
Zhongjin Lu
First Question
First Comment 
Can anything be done to solve the following problem? The dimension of x should be 3000. I have reduced it to 500 and Gurobi still cannot find the solution.


import gurobipy as gp
from gurobipy import GRB
import numpy as np
import scipy.sparse as sp
from scipy.sparse import coo_matrix
# Create a new model
m = gp.Model("matrix1")
# Set objective
obj = np.array(np.random.uniform(0, 3, 500))*0.01 # expected monthly returns
n = obj.shape[0]
# Create variables
x = m.addMVar(shape=n, lb=-GRB.INFINITY, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name="x")
lev = 4

# Create the auxiliary variable z to hold the result (z[i] = x[i] if x[i] < 0, else z[i] = 0)
zneg = m.addMVar(shape=n, lb=-GRB.INFINITY, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name="z")
# Create binary indicator variables for each x[i] < 0
b = m.addMVar(shape=n, vtype=GRB.BINARY, name="b")
# Indicator constraints for negative values
# Define the large constant M
M = 1000.0
for i in range(n):
    # If b[i] = 1, then x[i] < 0 -> z[i] = x[i]
    m.addGenConstrIndicator(b[i], True, zneg[i] == x[i], name=f"z_equals_x_if_negative_{i}")

    # If b[i] = 0, then x[i] >= 0 -> z[i] = 0
    m.addGenConstrIndicator(b[i], False, zneg[i] == 0, name=f"z_equals_0_if_positive_{i}")

    # Force b[i] to track whether x[i] is negative
    m.addConstr(x[i] <= M * (1 - b[i]), name=f"x_leq_constraint_{i}")  # b[i] = 1 if x[i] < 0
    m.addConstr(x[i] >= -M * b[i], name=f"x_geq_constraint_{i}")  # b[i] = 0 if x[i] >= 0


# Create an auxiliary variable to hold the absolute values of x
xabs = m.addMVar(shape=n, vtype=GRB.CONTINUOUS, name="z")
for i in range(x.shape[0]):
    # Add the absolute value constraint on the auxiliary variable
    m.addGenConstrAbs(xabs[i], x[i], name=f"xabs_constr_{i}")

# Add auxiliary variables to represent |x[i] - 1/7|
# Create the constant
constant_value = 1/n
# Create another auxiliary variable for x[i] - 1/7
aux = m.addMVar(shape=n, lb=-GRB.INFINITY, ub=GRB.INFINITY,  vtype=GRB.CONTINUOUS, name="aux")
zabs = m.addMVar(shape=n, vtype=GRB.CONTINUOUS, name="z")
for i in range(n):
    m.addConstr(aux[i] == x[i] - constant_value, name=f"aux_constr_{i}")
    # Add the absolute value constraint on the auxiliary variable
    m.addGenConstrAbs(zabs[i], aux[i], name=f"abs_constr_{i}")


m.setObjective(obj @ x - gp.quicksum(zabs)*0.25 + gp.quicksum(zneg)*2/12, GRB.MAXIMIZE) #50bps per trade; 50bps per month for shorting
# margin constraint
m.addConstr(gp.quicksum(xabs) <= lev , name="margin")
beta = np.random.uniform(0, 2, n)*5
Sigma_psd = np.diag(np.diag(np.random.uniform(4, 50, n))) + np.outer(beta, beta)

# Convert the dense Sigma_psd matrix into a sparse matrix using scipy
Sigma_sparse = coo_matrix(Sigma_psd)

# Define the upper bound V for the quadratic constraint;
V = 8 ** 2
# Add the quadratic constraint to the model
m.addConstr(x @ Sigma_sparse @ x == V, name="quadratic_constraint")

# Optimize model
m.optimize()
weights = x.X
print("Weights from Gurobi:", ", ".join(f"{w:g}" for w in weights))
print(f"Obj: {m.ObjVal:g}")


Related to

0

Comments

6 comments

Date Votes
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Hi  Zhongjin,

    Give our NoRel heuristic a go, by setting a non-zero value for NoRelHeurTime.  I think you'll be pleasantly surprised.

    I'm not sure I would hold out hope for getting a good bound (and therefore MIP gap) though, given it is a non-convex MIQCP.

    - Riley

    0
  • Zhongjin Lu
    First Question
    First Comment

    Riley,

    Thanks for the suggestion. If I want to get the closest approximation to the best possible solution within a certain amount of time, is there any other parameters I can set?

    Thanks,

    Zhongjin

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Hi Zhongjin,

    There are a lot of parameters but it is hard to guess what may work for your problem.  It sounds like you are less interested in the dual bound and more interested in good solutions, even if they are not provably good (via the dual bound).

    I would suggest running NoRel for a period of time, then switching to the main solve (root LP + tree) and setting MIPFocus=1, which tells the solver to focus on finding good solutions.  How long to run NoRel for depends on how well it is improving upon solutions, and how well solutions are improved in the tree.  When NoRel is successful it tends to produce solutions fast but the performance tapers off.  By monitoring your log you can get a feel for when the tapering happens.  One experiment you might want to try is run NoRel for 1,2,4,6,8 hours and save the solutions to MST files.  Then run experiments where these solutions are used as MIP starts with MIPFocus=1 (and NoRel) and see how well the solve does.

    - Riley

    0
  • Zhongjin Lu
    First Question
    First Comment

    Dear Riley,

    Following your suggestions, I tried to use some heuristic starting values. My question is that given that "Starting Obj Value using x[i].start" is 26.9437, why I got "User MIP start did not produce a new incumbent solution" when NoRel gives me substantially lower value: "Found heuristic solution: objective 4.7035734"

    How to guide Gurobi to search in bounds near the starting points? Thanks!

     

    # Create a new model
    m = gp.Model("matrix1")
    # Set objective
    obj = np.array(np.random.uniform(0, 3, 500)) # expected monthly returns
    n = obj.shape[0]
    # Create variables
    x = m.addMVar(shape=n, lb=-GRB.INFINITY, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name="x")
    lev = 4

    # Create the auxiliary variable z to hold the result (z[i] = x[i] if x[i] < 0, else z[i] = 0)
    zneg = m.addMVar(shape=n, lb=-GRB.INFINITY, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name="z")
    # Create binary indicator variables for each x[i] < 0
    b = m.addMVar(shape=n, vtype=GRB.BINARY, name="b")
    # Indicator constraints for negative values
    # Define the large constant M
    M = 1000.0
    for i in range(n):
        # If b[i] = 1, then x[i] < 0 -> z[i] = x[i]
        m.addGenConstrIndicator(b[i], True, zneg[i] == x[i], name=f"z_equals_x_if_negative_{i}")

        # If b[i] = 0, then x[i] >= 0 -> z[i] = 0
        m.addGenConstrIndicator(b[i], False, zneg[i] == 0, name=f"z_equals_0_if_positive_{i}")

        # Force b[i] to track whether x[i] is negative
        m.addConstr(x[i] <= M * (1 - b[i]), name=f"x_leq_constraint_{i}")  # b[i] = 1 if x[i] < 0
        m.addConstr(x[i] >= -M * b[i], name=f"x_geq_constraint_{i}")  # b[i] = 0 if x[i] >= 0


    # Create an auxiliary variable to hold the absolute values of x
    xabs = m.addMVar(shape=n, vtype=GRB.CONTINUOUS, name="z")
    for i in range(x.shape[0]):
        # Add the absolute value constraint on the auxiliary variable
        m.addGenConstrAbs(xabs[i], x[i], name=f"xabs_constr_{i}")

    # Add auxiliary variables to represent |x[i] - 1/7|
    # Create the constant
    constant_value = 1/n
    # Create another auxiliary variable for x[i] - 1/7
    aux = m.addMVar(shape=n, lb=-GRB.INFINITY, ub=GRB.INFINITY,  vtype=GRB.CONTINUOUS, name="aux")
    zabs = m.addMVar(shape=n, vtype=GRB.CONTINUOUS, name="z")
    for i in range(n):
        m.addConstr(aux[i] == x[i] - constant_value, name=f"aux_constr_{i}")
        # Add the absolute value constraint on the auxiliary variable
        m.addGenConstrAbs(zabs[i], aux[i], name=f"abs_constr_{i}")



    # margin constraint
    m.addConstr(gp.quicksum(xabs) <= lev , name="margin")
    beta = np.random.uniform(0.5, 2, n)*5
    Sigma_psd =  np.diag(np.random.uniform(10, 50, n)) + np.outer(beta, beta)

    # Convert the dense Sigma_psd matrix into a sparse matrix using scipy
    Sigma_sparse = coo_matrix(Sigma_psd)

    # Define the upper bound V for the quadratic constraint;
    V = 8 ** 2
    # Add the quadratic constraint to the model
    m.addConstr(x @ Sigma_sparse @ x == V, name="quadratic_constraint")
    m.setParam('NoRelHeurTime', 720)
    m.setParam('TimeLimit', 3600)     # Limit solving time to 10 minutes
    m.setParam('MIPGap', 0.01)       # Stop when a solution within 1% of the best bound is found
    # Set x equal to the vector "weights"
    Sigma_inv = np.linalg.inv(Sigma_psd)  # Inverse of Sigma
    weights = np.dot(Sigma_inv, obj)   # Dot product of inv(Sigma) and obj
    weights = np.sqrt(V/np.dot(weights.T, np.dot(Sigma_psd, weights)))*weights
    for i in range(len(weights)):
        x[i].start = weights[i]  # Set the starting value of x[i] to weights[i]

    # Print the objective value before optimization
    starting_obj_value = np.dot(weights.T, obj) - np.sum(np.abs(weights - constant_value)) * 0.25 + np.sum(weights[weights < 0]) * 2/12
    print(f"Starting Obj Value using x[i].start: {starting_obj_value:g}")

    # m.setParam('MIPFocus', 1)
    m.setObjective(obj @ x - gp.quicksum(zabs)*0.25 + gp.quicksum(zneg)*2/12, GRB.MAXIMIZE)  # 50bps per trade; 50bps per month for shorting

    # Optimize model
    m.optimize()

    # Save the NoRel heuristic solution to MST and write solution
    m.write("solution.mst")
    m.write("solution.sol")

    G_weights = x.X
    print("Weights from Gurobi:", ", ".join(f"{w:g}" for w in G_weights))
    print(f"Obj: {m.ObjVal:g}")

     

     

    0
  • Zhongjin Lu
    First Question
    First Comment

    Just one more piece of information: even the BestBd is below the obj value under the starting x values I give to Gurobi. This is puzzling to me. 

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Hi Zhongjin

    even the BestBd is below the obj value under the starting x values I give to Gurobi.

    This would imply the starting values are not feasible.

    If you fix the variables via lower and upper bounds the status is infeasible

    x[i].lb = weights[i]
    x[i].ub = weights[i]

    This article will contain some relevant advice:
    https://support.gurobi.com/hc/en-us/articles/16565310050833-Why-does-Gurobi-report-my-model-is-infeasible-when-it-has-a-feasible-solution

    - Riley

     

    0

Please sign in to leave a comment.

Didn't find what you were looking for?

New post