Skip to main content

Doubt in LP Constraints

Answered
V priya
  • Edited
Gurobi-versary
First Comment
First Question

I have the following code:

import numpy as np
import gurobipy as gp
from gurobipy import GRB

n = 3   
m = 2   
alpha = 0.9   

gamma = [0] * n
gamma[0] = 1
gamma = np.array(gamma)

Id = np.identity(n)

P_a1 = np.zeros((n, n))   
P_a2 = np.zeros((n, n))   

c = [60 + i * 5 for i in range(n - 3)]
c.extend([250, 120, 220])
c = np.array(c)

C_min = c + (alpha/(1-alpha))*min(c)
C_max = c + (alpha/(1-alpha))*max(c)

d_a1 = [65]*(n-3)
d_a1.extend([120, 200, 300])
d_a1 = np.array(d_a1)

d_a2 = [30 + i * 10 for i in range(n-2)]
d_a2.extend([150, 200])
d_a2 = np.array(d_a2)

D_min = np.minimum(d_a1, d_a2) + (alpha/(1-alpha))*min(np.minimum(d_a1, d_a2))
D_max = np.maximum(d_a1, d_a2) + (alpha/(1-alpha))*max(np.maximum(d_a1, d_a2))

ub = 50

for row in range(n-3):
    P_a1[row][row] = 0.2
    P_a1[row][row+1] = 0.8
P_a1[n-3][n-3] = 1
P_a1[n-2][0] = 0.8
P_a1[n-2][n-2] = 0.2
P_a1[n-1][n-1] = 1

for row in range(n-3):
    P_a2[row][row+1] = 0.3
    P_a2[row][n-2] = 0.6
    P_a2[row][n-1] = 0.1
P_a2[n-3][n-3] = 0.3
P_a2[n-3][n-2] = 0.6
P_a2[n-3][n-1] = 0.1
P_a2[n-2][n-2] = 1
P_a2[n-1][n-2] = 0.6
P_a2[n-1][n-1] = 0.4

U_lower_bar_a1 = np.zeros((n, n))   
U_upper_bar_a1 = np.zeros((n, n))   

U_lower_bar_a1[n-2][0] = -0.5
U_lower_bar_a1[n-2][n-2] = 0.3

U_upper_bar_a1[n-2][0] = -0.3
U_upper_bar_a1[n-2][n-2] = 0.5

U_lower_bar_a2 = np.zeros((n, n))  
U_upper_bar_a2 = np.zeros((n, n))  

for row in range(0, n-2):
    U_lower_bar_a2[row][0] = 0.5
    U_lower_bar_a2[row][n-2] = -0.5
U_lower_bar_a2[n-2][0] = 0.5
U_lower_bar_a2[n-2][n-2] = -0.7
U_lower_bar_a2[n-1][0] = 0.5
U_lower_bar_a2[n-1][n-2] = -0.3
U_lower_bar_a2[n-1][n-1] = -0.4

for row in range(0, n-2):
    U_upper_bar_a2[row][0] = 0.8
    U_upper_bar_a2[row][n-2] = -0.2
U_upper_bar_a2[n-2][0] = 0.7
U_upper_bar_a2[n-2][n-2] = -0.5
U_upper_bar_a2[n-1][0] = 0.7
U_upper_bar_a2[n-1][n-2] = -0.2
U_upper_bar_a2[n-1][n-1] = -0.3

M1 = np.vstack((np.identity(2*(n**2)), -np.identity(2*(n**2))))
M2 = np.hstack((np.hstack([np.identity(n)]*n), np.zeros((n, n**2))))
M3 = np.hstack((np.zeros((n, n**2)), np.hstack([np.identity(n)]*n)))
M4 = -np.vstack((M2, M3))

M_general = np.vstack((M1, M2, M3, M4))
M_polytope = np.vstack((M_general, np.random.rand(1, 2*(n**2))))

b_general = np.vstack(((U_upper_bar_a1.flatten('F')).reshape(-1, 1), (U_upper_bar_a2.flatten('F')).reshape(-1, 1), -(U_lower_bar_a1.flatten('F')).reshape(-1, 1), -(U_lower_bar_a2.flatten('F')).reshape(-1, 1), np.zeros((4*n, 1))))
b_polytope = np.vstack((b_general, np.random.rand(1)))
rho_dim = n*m

m = gp.Model('MR_problem_LP_approx_polytope')

rho = m.addMVar(rho_dim, lb=0.0, vtype=GRB.CONTINUOUS, name='rho')
z = m.addMVar(1, vtype=GRB.CONTINUOUS, name='z')
mu_c_pu = m.addMVar(len(b_polytope), lb=0.0, vtype=GRB.CONTINUOUS, name='mu_c_pu')
mu_d_pu = m.addMVar(len(b_polytope), lb=0.0, vtype=GRB.CONTINUOUS, name='mu_d_pu')

m.setObjective(z, sense=GRB.MINIMIZE)

m.addConstr( mu_c_pu.reshape(-1,1)@b_polytope <= z - np.tile(c,2)@rho + alpha*rho.reshape(-1,1)@np.vstack( U_lower_bar_a1, U_lower_bar_a2 )@(C_min - C_max).reshape(-1,1) )
m.addMConstr( M_polytope.transpose(), mu_c_pu, '=', alpha*np.append( rho[0:n]*n, rho[n, 2*n]*n )*np.tile( np.repeat( C_min, n ) ,2) )

In the last two constraints, I am getting error. I am new to gurobi in python. Please help !

Related to

0

Comments

5 comments

Date Votes
  • Gurobot (AI-generated response)
    Gurobi Staff Gurobi Staff

    The issue in your code seems to be related to the formulation of the constraints in Gurobi. Let's address them one by one.

    1. First Constraint Error: The first constraint appears to be incorrectly formed. Gurobi's addConstr method expects the left-hand side and right-hand side of the constraint to be linear expressions. In your code, however, the expression mu_c_pu.reshape(-1,1)@b_polytope seems to be a matrix multiplication, which might not be forming a valid Gurobi linear expression.

      To fix this, you need to ensure that the expression on each side of the <= operator is a linear expression. You can use Gurobi's quicksum or add linear expressions together using + and -.

    2. Second Constraint Error (addMConstr): The addMConstr method is used for adding multiple constraints at once. It seems like the dimensions of the matrices and vectors involved in the M_polytope.transpose(), mu_c_pu, '=' might not be matching appropriately.

      You need to ensure that the dimensions of M_polytope.transpose() and mu_c_pu are compatible for matrix multiplication and that the resulting dimension matches the dimension of the right-hand side of the constraint.

    To assist further, I would need to know the specific error messages you are receiving. This would provide clarity on whether it's a dimension mismatch, a type issue, or something else. In the meantime, here are general guidelines to troubleshoot:

    • Ensure that all the dimensions of matrices and vectors in your constraints match as required for matrix operations.
    • Make sure to use Gurobi's linear expression methods (LinExpr, quicksum, etc.) for forming constraints.
    • Check the types of all variables and expressions; Gurobi expects specific types for certain operations.

    If you can provide the exact error messages, I can offer more targeted advice.

    0
  • V priya
    Gurobi-versary
    First Comment
    First Question

    Regarding the first point, isn't it better to vectorize a constraint than use a for loop in terms of time ? 

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Hi,

    I think we've had a swing and a miss from our Gurobot here, and it's not aware that we can use addConstr with the Matrix API.

    When I run your code I'm told there is an error in this line:

    m.addConstr( mu_c_pu.reshape(-1,1)@b_polytope <= z - np.tile(c,2)@rho + alpha*rho.reshape(-1,1)@np.vstack( U_lower_bar_a1, U_lower_bar_a2 )@(C_min - C_max).reshape(-1,1) )


    An easy way of finding out where the problem is to add the following code before it, which simply executes parts of the code on separate lines:

    mu_c_pu.reshape(-1,1)@b_polytope
    np.tile(c,2)@rho
    alpha*rho.reshape(-1,1)
    alpha*rho.reshape(-1,1)@np.vstack( U_lower_bar_a1, U_lower_bar_a2 )
    alpha*rho.reshape(-1,1)@np.vstack( U_lower_bar_a1, U_lower_bar_a2 )@(C_min - C_max).reshape(-1,1)

    I'd expect the code to fail at one of these lines and it would give you a clue as to where the problem is, which can help you debug.

    In this case the first line fails.

    ValueError: Mismatching inner dimensions 1 vs. 49

    Now inspect the shapes of the objects being multiplied (either through a debugger, or interactively):

    >> mu_c_pu.reshape(-1,1).shape
    (49, 1)
    >> mu_c_pu.shape
    (49,)
    >> b_polytope.shape
    (49, 1)

    There's a few different "fixes" - including transposing, squeezing - but I'm wondering whether you need to reshape mu_c_pu at all.  If you're after a dot product then you can just use

    mu_c_pu@b_polytope

    Now this doesn't completely solve the issue - you have errors elsewhere in this expression but I'd encourage you to narrow down on the errors as shown above, and see if transposing or reshaping can resolve it.

    - Riley

    1
  • V priya
    Gurobi-versary
    First Comment
    First Question

    Hi 

    Thank you for the reply. I follow your help and see that:

    print( rho[0:n].shape) as 

    (3,) : means there are 3 columns (this is expected)
    and
    print( np.tile(rho[0:n],n).shape) as 
    (3,) : this is supposed to be 3*3=9. 
    Am i not allowed to use np.tile with decision variable of the model? But i define 
    it as m.addMVar, this means that numpy operations are applicable right ? 
    Also could you please provide answer to my previous question: 
     "Regarding the first point, isn't it better to vectorize a constraint than use a for loop in terms of time ?"
    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Hi,

    Our matrix API is built on top of numpy and scipy but this does not mean that numpy and scipy operations are compatible.  Sometimes they are compatible, but won't give the result you were hoping for.  Note, we have a few numpy-like functions in the gurobipy namespace: Matrix-Friendly API Functions, and of course there are the functions on the MVars themselves: 

    If you run the following code in an interactive Python session, such as Jupyter notebook:

    np.tile(rho[0:n],n)

    you will see that what you get back is a numpy array, where each element is a MVar - not exactly what you want.

    If you try and use gurobipy.hstack or gurobipy.concatenate on this result then you get back an MVar, but it has shape (9,).  If you try and specify axis=1 as an argument to gurobipy.concatenate then it will error since the MVars don't have this axis.  Since you want to stack it in the 2nd dimension you will need to explicitly reshape the elements to have this second dimension first.  Any of these will give you a 3x3 MVar:

    gp.concatenate(np.tile(rho[0:n].reshape(n,1),n), axis=1)
    gp.hstack(np.tile(rho[0:n].reshape(n,1),n))
    gp.concatenate([rho[0:n].reshape(n,1)]*n, axis=1)
    gp.hstack([rho[0:n].reshape(n,1)]*n)

     

    Also could you please provide answer to my previous question: 
     "Regarding the first point, isn't it better to vectorize a constraint than use a for loop in terms of time?"

    Yes, the vectorized approach should be faster than a for-loop, especially if matrices are sparse.

    - Riley

    0

Please sign in to leave a comment.

Didn't find what you were looking for?

New post