Skip to main content

MVar errors in Python

Answered
Taylor Petty
Gurobi-versary
First Question
First Comment

I have a fixed p and q and I want to minimize the sum of squared errors between A^T @ p and q with some restrictions. I.e., min_A (L2_norm(A^T p - q)), subject to A only having entries in {0,1,2} (preferably {0,1/2,1} but I couldn't figure out how to do that, so I figured out a workaround) and subject to the row sums of A, or column sums of A^T, each being exactly equal to 2.

I've tried it several different ways, but I can't get anything to work. The setup is the same regardless of my method:

p = np.array([0.6,0.1,0.3])
q = np.array([0.2,0.3,0.25,0.1,0.15])

m = gp.Model()

k = len(p)
J = len(q)

A = m.addVars((k,J),vtype=gp.GRB.INTEGER,lb=0,ub=2,name="A")

m.addConstrs((sum(A[j,:]) == 2 for j in range(A.shape[0])))

 

This is where it starts to diverge. I've tried following https://support.gurobi.com/hc/en-us/articles/360039343931-How-do-I-pointwise-multiply-a-numpy-vector-with-a-1-MVar-with-the-matrix-friendly-Python-API- and https://support.gurobi.com/hc/en-us/articles/360039587271-How-do-I-multiply-an-array-and-a-2-D-MVar-object-using-the-matrix-friendly-Python-API- and https://support.gurobi.com/hc/en-us/articles/360038943132 but I haven't figured it out yet.

 

One solution was to add auxiliary variables:

aux = m.addMVar(J)
quant = m.addMVar(1)
m.addConstrs((aux[j] == A[:,j]@p/2 - q[j] for j in range(J)))
m.addConstrs(quant == aux @ aux)

m.setObjective(norm,sense=gp.GRB.MINIMIZE)

But that third line of aux constraints returned GurobiError: Unsupported type (<class 'gurobipy.MLinExpr'>) for LinExpr addition argument.\

I also tried running that set of constraints in a loop like this:

for j in range(J):
    m.addConstr(aux[j] == A[:,j]@p/2-q[j])

But this returned the same error as before.

 

I also tried minimizing the sum of squares directly, and after some errors with pow() and ** I came up with this:

norm = sum((A[:,j] @ p/2 - q[j]) * (A[:,j] @ p/2 - q[j]) for j in range(J))

m.setObjective(norm,sense=gp.GRB.MINIMIZE)

But the norm definition returned ValueError: dimension mismatch, even though it's scalar multiplication. Using @ instead of * gave me GurobiError: Cannot multiply with an MLinExpr from the left.

Help, please!! I've been stuck for two days.

0

Comments

6 comments

Date Votes
  • Jaromił Najman
    Gurobi Staff Gurobi Staff
    • Edited

    You could try formulating your problem without the usage of MVars.

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

    m = gp.Model("test")
    p = np.array([0.6,0.1,0.3])
    q = np.array([0.2,0.3,0.25,0.1,0.15])
    K = len(q)
    J = len(p)

    A = m.addVars(K, J, vtype=GRB.INTEGER, lb=0, ub=2, name="A")
    # add constraints row_k(A) = 2
    m.addConstrs((gp.quicksum(A[k,j] for j in range(J)) == 2 for k in range(K)))
    # add constraints row_k(A)*p - q_k = aux_k
    aux = m.addVars(K, vtype = GRB.CONTINUOUS, name="aux")
    for k in range(K):
      m.addConstr(gp.quicksum(A[k,j]*p[j] for j in range(J)) - q[k] == aux[k], name="aux_constr_%d"%(k))
    # add constraint norm_aux = L2_norm(aux1,aux2,aux3)
    norm_aux = m.addVar(vtype=GRB.CONTINUOUS, name="norm_aux")
    m.addGenConstrNorm(norm_aux, aux, 2.0, "normconstr")
    m.setObjective(norm_aux, GRB.MINIMIZE)
    # write model to check whether it looks as expected
    m.write("myLP.lp")
    m.optimize()

    Note that I used the addGenConstrNorm method introduced with Gurobi v9.5.0 to model the \(L_2\) norm.

    Best regards,
    Jaromił

    1
  • Taylor Petty
    Gurobi-versary
    First Question
    First Comment

    Thank you so much for this! "addGenConstrNorm" was a great addition, too.

    If I let p vary in the problem, but leave everything else the same, it amounts to adding three lines:

    p = m.addVars(K,vtype=gp.GRB.CONTINUOUS,lb=0,ub=1,name="p")
    m.addConstr((gp.quicksum(p[i] for i in range(K))==1),name="psum")
    m.params.NonConvex = 2

    But then sometimes the model is infeasible (e.g. when q = [0.7, 0.2, 0.1], q = [0.79, 0.11, 0.1]) but sometimes it works just fine (q = [0.8, 0.1, 0.1] or [0.7, 0.3]). Is there something I'm missing about the solvability here? I'm not seeing why it would vary based on such small perturbations of q. (A slight note is that the norm uses A[i,j]/2 entries, and the columns of A are being dotted with p instead of the rows, e.g. m.addConstr(gp.quicksum(A[k,j]*p[k]/2 for k in range(K)) - q[j] == aux[j],name=f"aux_constr{j}").

    0
  • Jaromił Najman
    Gurobi Staff Gurobi Staff

    Did you try determining why the model is infeasible? See How do I determine why my model is infeasible?

    If this did not help, could you please post a minimal reproducible example where one model is feasible and the other is not?

    Best regards, 
    Jaromił

    0
  • Taylor Petty
    Gurobi-versary
    First Question
    First Comment

    I was unable to deduce anything meaninful from the m.computeIIS() output. However, I have realized that I'm only able to find a solution if an exact solution exists.

    The example below doesn't work. However, if p is set to [0.1,0.3,0.6] (or any permutation of those three numbers), it's feasible. This particular q and p came because I generated them using a known matrix and setting q = A^T p. This pattern held elsewhere -- if I set q and p to have an exact solution, it works. Perturbing them will result in failure.

    import gurobipy as gp
    import numpy as np
    q = np.array([0.1, 0.45, 0.3, 0.15, 0.])
    p = np.array([0.11,0.29,0.6]) #<-- change to [0.1,0.3,0.6] and it works
    K = len(p)
    J = len(q)
    m = gp.Model()

    A = m.addVars(K,J,vtype=gp.GRB.INTEGER,lb=0,ub=2,name="A")
    m.addConstrs((gp.quicksum(A[k,j] for j in range(J)) == 2 for k in range(K)))
    aux = m.addVars(J, vtype=gp.GRB.CONTINUOUS,name="aux")
    for j in range(J):
        m.addConstr(gp.quicksum(A[k,j]*p[k]/2 for k in range(K)) - q[j] == aux[j],
                  name=f"aux_constr{j}")
    norm_aux = m.addVar(vtype=gp.GRB.CONTINUOUS, name="norm_aux")
    m.addGenConstrNorm(norm_aux, aux, 2.0, name="normconstr")

    m.setObjective(norm_aux,sense=gp.GRB.MINIMIZE)
    m.optimize()

    However, of interesting note is if I use the same numbers (no exact solution) and I write out the norm expression explicitly, it seems to work, unless I've gotten something wrong (completely new to Gurobi here). Similarly to the first post using matrices, I'm not seeing how one expression is different from the other, or if one has advantages I'm misunderstanding.

    q = np.array([0.1 , 0.45, 0.3 , 0.15, 0. ])
    p = np.array([0.11,0.29,0.6]) 
    m = gp.Model()
    k = len(p)
    J = len(q)
    A = m.addVars(k,J,vtype=gp.GRB.INTEGER,lb=0,ub=2,name="A")
    for row in range(k):
        m.addConstr((sum(A[row,j] for j in range(J))==2),
                     name=f"row{row}sum")
    norm = gp.quicksum(pow(gp.quicksum(A[i,j]*p[i]/2 for i in range(k))-q[j],2)
                       for j in range(J))
    m.setObjective(norm,sense=gp.GRB.MINIMIZE)
    m.optimize()
    0
  • Jaromił Najman
    Gurobi Staff Gurobi Staff

    Note that variable lower bounds are \(0\) per default. However, the values of your \(\texttt{aux[i]}\) variables can be negative. Allowing for negative lower bounds

    aux = m.addVars(J, lb=-gp.GRB.INFINITY, vtype=gp.GRB.CONTINUOUS,name="aux")

    solves the issue.

    Best regards, 
    Jaromił

    1
  • Taylor Petty
    Gurobi-versary
    First Question
    First Comment

    That was it! Thank you so much!

    0

Please sign in to leave a comment.

Didn't find what you were looking for?

New post