Issue in indexing gurobi matrix variable in constraints in python
回答済み
Hi,
I am working on an optimization problem. I have an issue in indexing the gurobi matrix variable in the constraint part. I am unable to solve the matrix dimension issue. Could anybody have a look at my script and guide me how to fix it?
#In the nested loop of the constraint part:
Ts*np.linalg.inv(M)@Bu@dDelta)@useq[i-k-2]+s1[:,k]
I am getting incompatible dimension error because of s1 variable. However, I am slicing the s1 matrix 3by1 which is of size 3x5 already defined. I don't understand where I am wrong.
Thank you very much in advance.
#python code
N=5
m=gp.Model("matrix")
Ts=0.2
M=np.eye(3)
D=np.eye(3)
C=np.eye(3)
Bu=np.array((3,4)
pos=np.array([0,0,0]).T
u=np.array([0,0,0,0]).T
xseq=[]
xseq.append(pos)
useq=[]
#1st loop
for i in range (1,N+1):
xseq.append(pos)
useq.append(u)
#2nd loop
for i in range(N+1,7):
s1=m.addMVar((3,N)) #s1 matrix of 3x5
A=np.eye(3)
deltaBut1=m.addVar(1)
deltaBut2=m.addVar(1)
deltaBut3=m.addVar(1)
deltaBut4=m.addVar(1)
deltaBut=np.array([deltaBut1,deltaBut2,deltaBut3,deltaBut4]).T
dDelta=np.diag(deltaBut)
obj1=0
for k in range(0,N): # N=5
m.addConstr(xseq[i-k-1]-A@xseq[i-k-2]==(Ts*np.linalg.inv(M)@Bu@dDelta)@useq[i-k-2]+s1[:,k])
obj1=obj1+s1[:,k]@s1[:,k]
m.setObjective(10**3*(obj1))
m.optimize
-
I would like to add one more thing that:
m.addConstr(xseq[i-k-1]-A@xseq[i-k-2]==(Ts*np.linalg.inv(M)@Bu@dDelta)@useq[i-k-2]+s1[:,k])
the left side gives 3x1 shape while right side gives MLinExpr of shape (3,). Could anybody help me?
0 -
Hi,
The \(\texttt{Bu}\) object is only of dimension 1x2. Additionally, it might be easier to work with \(\texttt{s1}\) as Var instead of an MVar object (see your other post).
Best regards,
Jaromił0 -
Thanks Jaromil. Sorry it was my mistake. Actually Bu object is 3x4 matrix. Using s1 as AddVars the problem solved in setting the objective function. However, I am still having difficulty in setting the constraints.
s1=m.addVars(3,N)
for k in range(0,N): # N=5
for j in range(0,3):
obj1=obj1+s1[j,k]*s1[j,k] #setting objective
#adding constraints at each iteration of k i.e 0,1,2,3,4
m.addConstr(xseq[i-k-1]-A@xseq[i-k-2]==(Ts*np.linalg.inv(M)@Bu@dDelta)@useq[i-k-2]+s1[:,k])Jaromił Najman Could you please guide me on how to modify the constraints part to get the s1 all rows and kth column at each iteration of k? Thank you very much.
0 -
I am not sure what you want to express with this constraint.
The left hand side of the constraint is a 1x3 double array. The part
Ts*np.linalg.inv(M)@Bu@dDelta)@useq[i-k-2]provides a 1x3 array with linear expressions and then you are trying to add a 1x3 variable array. So what you are trying to express is
[1 2 3] = [x y z] + [s1_1 s1_2 s1_3]
So do you actually want to write 3 constraints? Or do you want to sum everything and actually express
1+2+3 = x+y+z + s1_1+s1_2+s1_3
Best regards,
Jaromił0 -
I want to get this.
[1 2 3] = [x y z] + [s1_1 s1_2 s1_3]
I need three constraints not their summation.
I need to achieve :
m.addConstr([1 2 3]-[x y z]==[s1_1 s1_2 s1_3]) #at each iteration of k loop
0 -
I understand. What you could go for is
for k in range(0,N): # N=5
for j in range(0,3):
obj1=obj1+s1[j,k]*s1[j,k]
# both terms have dimension 3
term1 = xseq[i-k-1]-A@xseq[i-k-2]
term2 = (Ts*np.linalg.inv(M)@Bu@dDelta)@useq[i-k-2]
for l in range(0,3):
m.addConstr(term1[l] - term2[l] + s1[l,k] == 0)Best regards,
Jaromił0 -
Jaromił Najman Thank you very much for your kind help. Yes this is what I wanted. However, this constraint gives me "model is infeasible or unbounded" error.
0 -
You could use the write() function to write an LP file which you can then analyze and check whether the problem actually looks like what you expect it to be
m.write("myLP.lp")With the InfUnbdInfo parameter, you can tell Gurobi to try to determine whether the model is unbounded or infeasible.
The Knowledge Base article How do I determine why my model is infeasible? might also be helpful.
Best regards,
Jaromił0 -
Jaromił Najman Thanks a lot. By rearranging the constraints as
m.addConstr(term1[l] == term2[l] + s1[l,k])
I get the correct constraints i.e 0.
Now its working.
0 -
Yes, correct, I made a mistake when writing down the constraint. Good to hear that it works now!
Best regards,
Jaromił0 -
Hi,
I would like to ask one question. My algorithm has two optimization objectives. In the first loop, I do optimization of the objective 1 and save the records. In the next loop, first I add constraints to the objective 2 based on the N values. Then I solve the objective 2 and keep the record. At the same iteration I need to solve the objective 1 that uses delta optimized value from the objective 2 optimization. I described the scenario in the following algorithm.
My question is that would it be possible to use the same model i.e m=gp.Model("test") for both optimization problem. Currently I tried to solve like this but I am not sure its correct or not. When I run optimization of one objective at a time, the model is ok but when I run the code as I described in the algorithm, gurobi model stops working.
I am very beginner in the programming. Could you please check the algorithm and suggest the right way to follow? I really appreciate your kind responses. Thank you very much.
for i in range(1,N+1):
#do optimization of objective 1
for i in range(N+1,700):
for k in range (0,N):
#add constraints to objective 2
#do optimization of objective 2
#keep record of objective 2
.
.
#do optimization of objective 1 #Note: here objective 1 uses delta from objective 2
#solve the system
0 -
Hi,
The idea you described should work.
When exactly does Gurobi stop working?
You can work with the same model. You have to make sure that you set the correct objective via the setObjective function.
Then, once you have optimized for the first objective, added new constraints and solved for the second objective, you can remove old constraints using the remove function. It is best to keep track of all constraints you want to remove. In you case you could have a list of constraints, which you add and then remove then before solving for objective 1.
for i in range(N+1,700):
constrs = []
for k in range (0,N):
constrs.append(m.addConstr(...))
m.setObjective(obj2)
# optimize objective 2
# keep record of objective 2
# remove constraints used for objective 2
m.remove(constrs)
#do optimization of objective 1 #Note: here objective 1 uses delta from objective 2
m.setObjective(obj1)
#solve the systemIs this what you were looking for?
Best regards,
Jaromił0 -
Hi,
Thank you very much for your help. I was doing this. Unfortunately, the gurobi model does not work with objective 2. When I run objective 1 both in first and second loop. It is working properly. However, when I integrate the objective 2 in the loop. The model shows infeasible solution. I am reviewing the whole script. Hope I will get error. Thanks again.
0 -
You might be interested in the Knowledge Base article How do I determine why my model is infeasible?
Best regards,
Jaromił0 -
Hi,
When I was trying to get the gurobi var after running objective 2. The model was showing infeasible solution and no attribute with X.
Surprisingly when I set condition:
if m.status==GRB.OPTIMAL:
deltabut1.XI was able to get the gurobi var values and the model was working perfectly. Another thing I added the variable names during declaration. Now I am getting optimal solutions for both objective functions.
I would like to thank you again and greatly appreciated your kind responses.
0 -
Hi,
Great to hear that you models work now. Please note that the code
if m.status==GRB.OPTIMAL:
deltabut1.Xworks, because you can only access the solution values X when at least one feasible solution has been found. This does not mean that all your models are feasible. You could check for infeasible models via
if m.status==GRB.INFEASIBLE:
print("model is infeasible!")
if m.status==GRB.INF_OR_UNBD:
print("model is infeasible or unbounded!")However, if you are only interested in the models which are solved to optimality and you are free to ignore all other solution outcomes, then your approach is correct.
Best regards,
Jaromił0 -
yes you are right.
When I catch the infeasible or unbounded status, the model returns the following "model.ilp" after computing m.computeIIS(). This means that after some iteration the model is infeasible with the following constraints
Subject To
con3: d3 <= 1
con4: d4 <= 1
con9: C0 + 3.90347e-06 d1 + 3.73644e-06 d2 - 2.95246e-06 d3
- 2.84865e-06 d4 = 3.27893e-11
Bounds
d1 >= 1
d2 >= 1
d3 free
d4 free
Constant = 1
EndHowever, I need this from the model.
In the script the delta variables (d1,d2,d3,d4) are subject as : 1>=deltaBut>=0
m.addConstr(deltaBut1<=1,name='con1')......m.addCosntr(deltaBut4>=0,name='con8')
#and also the "con9" should be 0 according to the constraints
m.addConstr(term1[l] == term2[l] + s1[l,k])do you have any idea on this problem? why the model is not following the defined constraints?
Kind regards,
wa
0 -
Hi wa,
You state that all your \(d_i\) variable have to between in \([0,1]\). In your previous comment, you define your \(d_i\) variables via
deltaBut1=m.addVar(1)
[...]This sets the lower bound of the variable to 1. In order to set the upper bound to 1, you have to provide the argument's name, i.e.,
deltaBut1=m.addVar(ub=1)
[...]Note that the default lower bound for optimization variables is 0.
If this does not help, you can use the feasRelax function to construct a feasibility relaxation and determine the change you have to perform in order to make your problem feasible.
Best regards,
Jaromił0 -
I hope you are doing very well. I would like to ask one question.
In python 2.7, the @ operator is not supported.
What is the alternative of this if I want to set my objective as I did with python 3.7. With python 2.7 is gives error with (u*R*u)
Thanks a lot!
m.setObjective(u@R@u)
0 -
Hi wa,
The best solution would be to update to a higher Python version. However, this may not always be possible due to different reasons. In this case, you have to compute the matrix multiplications by hand, i.e., with \(\texttt{for}\)-loops, using quadratic expressions. To easily access the entries of variables \(u\), you could introduce the variables as normal variables and not MVars or use the tolist function to generate a list of your MVars and then work with this variable list.
Best regards,
Jaromił0 -
thank you very much. In my case I have to adopt Python 2.7. I used tolist function and it worked in setting objective function.
u=m.AddMvar(shape=4)
uu=u.tolist()
obj1=np.matmul(uu,R)
objective=np.matmul(obj1,uu)
m.setObjective(objective)But I also have to set constraints in python 2.7 version format.
m.addConstr(Ts*np.linalg.inv(M)@Bv@tau==Ts*np.linalg.inv(M)@But@u)
In the constraints part, the left side I get easily by just simply multiplying the matrices. I need to achieve from right side:
right_con1=np.matmul(Ts*np.linalg.inv(M),But)
right_con2=np.matmul(right_con1,uu)However here, the u.tolist() is not working?
Could you please check?
Thank you a lot.
Kind regards,
WA
0 -
Hi WA,
Could you please post a short snippet of how you are using the tolist function and what exact error you get?
Best regards,
Jaromił0 -
Here is my script.
This generates infeasible model because of uu variable multiplication in right_term12.
With the following the model is Ok. But I need to get rid of @ operator in python 2.7
m.addConstr(Ts*np.linalg.inv(M)@Bv@tau==Ts*np.linalg.inv(M)@But@u)
Kind regards,
WA
0 -
Hi WA,
Since you have to avoid the @ operator anyway, you could work with Vars instead of MVars. This way, you would avoid having to work with the tolist function.
Could you post a whole reproducible code snippet? In your screenshot, there are terms missing, \(\texttt{Ts, M, Bv, old\_tau, But}\).
The safest way to avoid the @ operator is working with \(\texttt{for}\)-loops, i.e., implementing the matrix multiplication yourself. This means more work but will work among all Python versions including 2.7.
Best regards,
Jaromił0 -
Hi Jaromil,
I would like to thank you for your kind help and support for all my queries.
I did not bother you again because I found a solution for Python version and I was able to run my script with Python 3.
Your help is highly appreciated.
Kind regards,
WA
0
サインインしてコメントを残してください。
コメント
25件のコメント