Quadratic Expression

Hi Everyone,

I am trying to solve the following code, but it takes forever the solve the model even tough it is not a big problem. Can you help me what is the problem? I believe I cannot handle the quadratic expression parts. I have l2 norm in both objective function and constraints.

Thank you!

import numpy as np
from math import *
from gurobi import *
from scipy.stats import *

A=np.array([[1,1,0,0,1,1],
            [0,0,1,1,1,1]])
price=np.array([150, 100, 120, 80, 250, 170])
c=np.array([90,90])
# expected_demand=np.array([30,60,20,80,30,40])
expected_demand=np.array([30,60,20,60,30,40])
resources=len(A)
product=len(A[0])
P=np.identity(product)
SD=np.array([5.47,7.74, 4.47,8.94,5.47,6.32])
# SD=np.array([5.47,20, 4.47,8.94,5.47,13])
for i in range(product):
    for j in range(product):
        if i==j:
            P[i,i]=SD[i]

epsilon=1    #epsilon gamma
b=0.3        #rho
omega=3      #l2 norm omega

resources=len(A)
product=len(A[0])

M=Model("box-ellipsoid Uncertainity Network Optimization")

y= M.addMVar(product, vtype = GRB.CONTINUOUS, lb = 0, ub = 1, name = 'y')
x= M.addMVar(product, vtype = GRB.CONTINUOUS, name = 'x')
sqrt_expr1 = M.addMVar(1, vtype=GRB.CONTINUOUS, name='sqrt_expr1')
sqrt_expr2= M.addMVar(resources, vtype=GRB.CONTINUOUS, name='sqrt_expr1')
x_const={}
for i in range(resources):
    for j in range(product):
        x_const[i,j]=M.addVar( vtype = GRB.CONTINUOUS, name = f'obj_var_dual_{i,j}')
pi_constr_var1=M.addMVar(resources, vtype = GRB.CONTINUOUS, lb = 0, name = 'pi_constr_var1')
obj_var1=M.addVar( vtype = GRB.CONTINUOUS, lb=0, name = 'obj_var1')
tau=M.addMVar((resources,product), vtype = GRB.CONTINUOUS, lb = 0, name = 'tau')
w=M.addMVar(product, vtype = GRB.CONTINUOUS, lb = 0, name = 'omega')

#objective function
expr1 = quicksum(x[j] * x[j] for j in range(product))
M.addConstr(sqrt_expr1 * sqrt_expr1 == expr1)
objFun = quicksum(price[j] * y[j] * expected_demand[j] for j in range(product)) - omega * sqrt_expr1 - epsilon * obj_var1
M.addConstr(sum(w[j] for j in range(product))<=obj_var1)
M.addConstrs(w[j] >= (price[j] * y[j] * P[j,j])-x[j] for j in range(product))
M.addConstrs(w[j] >= (-price[j] * y[j] * P[j,j])+x[j] for j in range(product))
M.update()
#constraints
expr2=[0]*resources
for i in range(resources):
    expr2[i] = quicksum(x_const[i,j]*x_const[i,j] for j in range(product))
    M.addConstr(sqrt_expr2[i] * sqrt_expr2[i] == expr2[i])
    M.addConstr(c[i] >= quicksum(A[i][k] * y[k] * expected_demand[k] for k in range(product))+omega*sqrt_expr2[i]+epsilon*pi_constr_var1[i])
    M.addConstr(pi_constr_var1[i] >= sum(tau[i,j] for j in range(product)))
M.addConstrs(tau[i,j] >= P[j,j]*y[j]*A[i,j]-x_const[i,j] for i in range(resources) for j in range(product))
M.addConstrs(tau[i,j] >= -(P[j,j]*y[j]*A[i,j]-x_const[i,j]) for i in range(resources) for j in range(product))
M.update()
M.setObjective(objFun, GRB.MAXIMIZE)
M.setParam('OutputFlag', False)
M.setParam('NonConvex', 2)

M.update()
M.optimize()
status = M.status

values = y.X
M.update()

y1=values
print(y1)

Jaromił Najman

Gurobi Staff

January 16, 2024 08:26 Edited

Hi Irem,

Could you please try using the addGenConstrNorm method to introduce the L2 norm?

Are you using the latest Gurobi version 11.0.0? If no, please try updating to the latest version if possible.

Could you please also share a log snippet showing the optimization process?

I just tried your code on my personal Mac with Gurobi v11 and it reaches a MIPGap of ~0.7% in only a few seconds. What is the required solution quality for your model, i.e., what is a satisfying MIPGap value for your application?

I can also see the warning

Warning for adding constraints: zero or small (< 1e-13) coefficients, ignored

This means that some of the coefficients you compute are very small. If these coefficients are of relevance, you would have to re-scale your model. If they are not of any relevance, I would recommend removing them from the start.

You should also try adding upper bound for your variables. I don't know the details of your model, but I would guess that an upper bound of 1 might be fitting for most of your variables.

Best regards,
Jaromił

İrem Bahtiyar

January 16, 2024 09:09

Mr. Najman Hi,

Thank you for your response, I appreciate it! I did not specify MIPGap at first, but now I did and it solved my problem quickly! MIPGap of 0.03 is enough for me. Thank you!

I am aware of that warning but I could not understand how can I get rid of it. I would appreciate elaboration on this.

Best,

İrem

January 16, 2024 09:28

Hi İrem,

It may not be easy to detect. One way to find out would be to print every coefficient you compute just to find out which ones are so tiny.

For example, you could print the values of \(\texttt{price * P}\)

for j in range(product):
  if price[j] * P[j,j] < 1e-6:
    print(f"index {j} value: {price[j] * P[j,j]}")

You would have to do this manually for other coefficients you generate. This way, you can detect which coefficients are very small and why and possibly remove them from the model completely. It is usually OK to remove terms with tiny coefficients. However, please note that you can completely remove these terms from your model only if the tiny coefficient don't play a significant role on the model and solution behavior. If these rows are required for the correctness of the application, you should have a look at our Guidelines for Numerical Issues and try to re-scale the model.

January 16, 2024 09:31

I understand, let me check my model. Thank you for your help!

Quadratic Expression

Comments

Didn't find what you were looking for?