No sensible solution but Code dont have errors
Answered
Hi, I wrote a code in Gurobi Python API for minimizing production costs which have a variable and a fixed cost part. There are a number of tasks which are assigned to a number of stations whereby cycle time should not be exceeded. Also 4 types of machines are given which the program should evaluate which ones are being used, the machines with higher var. costs have higher fixed costs and vice versa. The part with different types of machines is making me the most headaches, because the duration and the cost of assigning a task is not dependent on the station indices but only on the machine typ. I dont know the problem but the program gives no errors but also not really a solution is calculated. I think the problem might be that the constraints are not programmed correctly so that they dont have a impact on the Optimization Model. Also it could be that my objective function is not on target. Following is my code with reduced tasks, stations etc.:
import gurobipy as gp
from gurobipy import GRB, quicksum
stations = [1,2,3,4]
tasks = ['w01','w02','w03']
durations = {('w01','BTC1'):200,
('w02','BTC1'):500,
('w03','BTC1'):300,
('w01','BTC2'):100,
('w02','BTC2'):250,
('w03','BTC2'):150 }
successors = {'w01':[],'w02':['w01'],'w03':['w01']}
types = ['BTC1','BTC2']
arcs = [(i,j) for i in tasks for j in types]
arcs2 = [(i,j) for i in stations for j in tasks]
C = 500
C_min = 100
edges, v_cost, f_cost = gp.multidict({
('BTC1','w01'): [150,200000],
('BTC1','w02'): [200,200000],
('BTC1','w03'): [180,200000],
('BTC2','w01'): [120,300000],
('BTC2','w02'): [160,300000],
('BTC2','w03'): [100,300000]
})
m = gp.Model()
y = m.addVars(stations, vtype=GRB.BINARY, name='y')
x = m.addVars(arcs, vtype=GRB.CONTINUOUS, name='x',lb=0,ub=1) # Variable for share of assignment between 0 and 1
a = m.addVars(arcs2, vtype=GRB.CONTINUOUS, name='a',lb=0,ub=1)
# every task is exactly done by a share of 1, so every task is fully completed over the number of stations
m.addConstrs(quicksum(a[i,j] for j in tasks) == 1 for i in stations)
#Capacity constraint
m.addConstrs(quicksum(durations[i,j]*x[i,j] for i in tasks for j in types)
<= C*y[j] for j in stations)
# there is also a minimum Capacity
m.addConstrs(quicksum(durations[i,j]*x[i,j] for i in tasks for j in types)
>= C_min*y[j] for j in stations)
# precedence constraints
m.addConstrs(quicksum(p*x[k,j] for j in types for p in stations) >=
quicksum(p*x[i,j] for j in types for p in stations)
for k in tasks for i in successors[k] if k!= 'w01')
m.setObjective(x.prod(v_cost)+y.prod(f_cost), GRB.MINIMIZE)
m.optimize()
m.optimize() only returns:
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 14 rows, 22 columns and 76 nonzeros
Model fingerprint: 0x4cbb069e
Variable types: 18 continuous, 4 integer (4 binary)
Coefficient statistics:
Matrix range [1e+00, 5e+02]
Objective range [0e+00, 0e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Presolve removed 11 rows and 16 columns
Presolve time: 0.01s
Presolved: 3 rows, 6 columns, 14 nonzeros
Variable types: 6 continuous, 0 integer (0 binary)
Root relaxation: objective 0.000000e+00, 0 iterations, 0.01 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
* 0 0 0 0.0000000 0.00000 0.00% - 0s
Explored 1 nodes (0 simplex iterations) in 0.12 seconds (0.00 work units)
Thread count was 4 (of 4 available processors)
Solution count 1: 0
Optimal solution found (tolerance 1.00e-04)
Best objective 0.000000000000e+00, best bound 0.000000000000e+00, gap 0.0000%
It says optimal solution found but obviously there is not really a solution. I would really appreciate your help.
-
How exactly can I insert the solution point values into my constraints and check if they are right and doing the things they should? I am a beginner so I would really appreciate your further help.
Your model is small enough such that you can do this easily by hand. Note that all your constraints will be fulfilled because your model is feasible. The question you should ask yourself now is, are those solution point values the ones you would expect. For example, are the stations \(y\) open that you would expect to be open in an optimal solution? What are the arc values? Do they have the values, they should have in an optimal solution?
Do you have a model formulation written on paper somewhere? You could compare the LP file with the formulation it should depict and check whether every constraints looks exactly as it should.
1 -
To debug such issues, it is best to write the solution file and the model file after optimization has finished.
m.optimize()
m.write("solution.sol")
m.write("mymodel.lp")and then inspect it. This means, that you should check whether the solution point values make sense. You can then insert the solution point values into your constraints and see which constraints are not doing what they actually should do. The \(\texttt{mymodel.lp}\) file allows you to analyze whether the model is formulated correctly.
Best regards,
Jaromił0 -
Thanks for your help Sir. I did export the LP file and the .sol file but I am not sure about the next step you described.: "You can then insert the solution point values into your constraints and see which constraints are not doing what they actually should do"
How exactly can I insert the solution point values into my constraints and check if they are right and doing the things they should? I am a beginner so I would really appreciate your further help.
0 -
It's odd that the coefficient statistics show
Objective range [0e+00, 0e+00]
This suggests that all of the objective coefficients are zero, which would explain why an optimal objective value of 0.000000000000e+00 is found.
0 -
Okay thats a good hint Robert Fourer. Maybe I restricted the objective somehow but did not notice how. Is there any suspicious line within the code where you could think could cause some kind of restriction? Of course I dont want the objective to be in a interval of [0,0]. Thanks in advance
0 -
Jaromił Najman I did read the LP of my code, huge thanks for that hint because all of the constraints where really of place and did not do what I wanted mathematically. The one Constraint Unfortunetely I could not manage to get right or even understand the LP file to be honest. What I want are precedence constraints and the precedences should be like that:
{'w01': [],'w02': [],'w03': [],'w04': [],'w05': ['w03'],'w06': ['w01'],
'w07': ['w02'],'w08': ['w04'],'w11': ['w05'],'w09': ['w04'],'w10': ['w06'], 'w12': ['w07'],'w14': ['w08'],'w13': ['w09'],'w15': ['w10'],'w16': ['w11'], 'w17': ['w13'],'w18': ['w15'],'w19': ['w16'],'w20': ['w17']}What I provided is the Output of my successors dictionary.' Please note that it is read like task': [predecessor]
The mathematic model of the constraints is that:
But one addition is that I not only have a number of stations but also different types of machines that are assigned to the stations and deviate in form of costs and durations. To implement this in Gurobi I did use this code:
m.addConstrs(quicksum(p*x[k,j] for j in types for p in stations) >=
quicksum(p*x[i,j] for j in types for p in stations)
for k in tasks for i in successors[k] if k!= 'w01')In the LP of my code the constraints are listed like that (dont wonder I added 2 types of machines from original code but this should not change the problem):
R28: - 15 x[w03,BTC1] - 15 x[w03,BTC2] - 15 x[w03,BOC1] - 15 x[w03,BOC2]
+ 15 x[w05,BTC1] + 15 x[w05,BTC2] + 15 x[w05,BOC1] + 15 x[w05,BOC2]
>= 0and so on. Is the constraint that I want to implement right in comparision to the R28 line of the LP?
Thanks in advance and sorry for long text
0 -
To better recognize your constraints, you should always give them a name, e.g.,
m.addConstrs(quicksum(p*x[k,j] for j in types for p in stations) >=
quicksum(p*x[i,j] for j in types for p in stations)
for k in tasks for i in successors[k] if k!= 'w01', name="myconstr")This should help you in fixing the constraint. If you don't manage to fix it, could you please provide a smaller example concentrating only on the one constraint you are trying to implement? This would make tracking down the error much easier.
Regarding the objective function. The dictionaries \(\texttt{v_cost}\) and \(\texttt{f_cost}\) do not fit the \(x,y\) variables.
Your \(x\) variables are defined over arcs which is given as \(\texttt{tasks x types}\) but the \(\texttt{v_cost}\) dictionary is defined over \(\texttt{types x tasks}\). The key value has to be exactly the same in both dictionaries, cf. tupledict.prod() documentation. To fix the \(x\) part, you can write
edges, v_cost, f_cost = gp.multidict({
('w01','BTC1'): [150,200000],
('w02','BTC1'): [200,200000],
('w03','BTC1'): [180,200000],
('w01','BTC2'): [120,300000],
('w02','BTC2'): [160,300000],
('w03','BTC2'): [100,300000]
})For \(y\), you will have to define a different \(\texttt{f_cost}\) list because your \(y\) variables are defined over stations \(\{1,2,3,4\}\) while the current \(\texttt{f_cost}\) list is defined over the 6 edges
('BTC1','w01'), ('BTC1','w02'), ('BTC1','w03'), ('BTC2','w01'),
('BTC2','w02'), ('BTC2','w03')One possible solution would be
f_cost = {1: 200000, 2: 200000, 3:300000, 4:300000}but I don't know which values should be used here.
0 -
Jaromił Najman thanks for your active help. I managed to get all of the constraints like I want them to be as I looked up at the LP file. Now for the cost the issue is that they are not dependent on stations. I should explain the problem more specific. There are different types of machines from which a company can choose. Each has advantages and disadvantages, to brief it up shortly the machines with high investment costs do have lower variable costs. I used the stations array only as a supporting variable. When only define the 2 types of machines the programm assumes that there are only 2 machines which can be assigned to tasks. But there are not only 2 machines, there are only 2 types of machines from which the solution can choose which ones are used how often. E.g. it could be possible that machine type 1 is used 3 times and machine type 2 4 times, but I cannot do that because if I delete the stations array the program thinks there are only 2 machines available. Therefore I introduced the stations which have no impact to the costs. Therefore they are not relevant and theoretically not needed but I dont know how to deal without them. Any idea on this topic? Thanks
0 -
Okay I fixed the x variable but with the order like the edges e.g. ('BTC1','w01') just because the notation in mathematical model is also like that. But I dont think that it is making any difference. Also i get know a Objective range in the coefficients. The solver returns now:
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 44 rows, 184 columns and 396 nonzeros
Model fingerprint: 0x56604daf
Variable types: 180 continuous, 4 integer (4 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+03]
Objective range [2e-02, 3e+05]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Presolve removed 20 rows and 100 columns
Presolve time: 0.01s
Presolved: 24 rows, 84 columns, 296 nonzeros
Variable types: 80 continuous, 4 integer (4 binary)
Found heuristic solution: objective 0.0000000
Explored 0 nodes (0 simplex iterations) in 0.02 seconds (0.00 work units)
Thread count was 4 (of 4 available processors)
Solution count 1: 0
Optimal solution found (tolerance 1.00e-04)
Best objective 0.000000000000e+00, best bound 0.000000000000e+00, gap 0.0000%Also I should mention that I rewrote the y variable to be dependent on types and not on stations anymore like that:
y = m.addVars(types, vtype=GRB.BINARY, name='y')
Another news is that my Objective Function is now shown in the LP which was not the case before. It looks like that:
Minimize
200000 y[BTC1] + 300000 y[BTC2] + 150000 y[BOC1] + 180000 y[BOC2]
+ 0.8964 x[BTC1,w01] + 0.747 x[BTC1,w02] + ... +The values for variable as well for fixed costs are all as they should be.
But still I get no result which is really annoying. I think there is really something with the objective function wrong that it doesnt count the costs properly. I will be very happy for further help because I am near to the solution now. Thanks
0 -
But I dont think that it is making any difference.
It is not making a significance on the formulation but it does affect the programming part. As you noticed, now the \(x\) variables are present in the objective function.
Did you have a look at the solution point? Since the objective value is \(0\), all \(x,y\) variables have to be \(0\) in the current optimal solution. Probably, some of your constraints have a bug and are currently not binding.
0 -
So in the solution point I recognized something that I think might bring the error. I defined a decision variable a with arcs2 values given to it. And in the solution some of the a´s are 1 so they have values but I dont want them to get values I want the x to get values. Therefore the x and y variables as you sad are all equal to 0.
The main problem for me is: I still dont know how I can formulate this without the stations array because it disrupts me in finding a solution. I think I could just make a 3-tuple. E.g. instead of having ('BTC1','w01') I could say ('BTC1','w01',1) 1 for station 1 and so on. But it is not really an efficient solution I think.
So that you really understand my problem in the right way:
I just have the stations variable to say that there can be stations in different types of machines, e.g. BTC1 can be used twice and BOC1 never, because that is what the optimal solution finds for the minimum costs. I dont want to be restricted to must having each one of the machine types. Is there a way that I can do this without the help of stations and therefore also without the a variable. In this case I think the values for a would come over to the values of x and there will be a solution.
Thanks so much
0 -
Jaromił Najman so I think i fixed all of the problems now and when I look at the LP everything is just fine.
I did define the y variable as tuple of (' station ',' type ') and it solves my problems. In the solution point I get that: y[1,BTC1] 1y[1,BTC1] 1
y[1,BTC2] 1
y[1,BOC1] 1
y[1,BOC2] 1
y[2,BTC1] 1
y[2,BTC2] 1
y[2,BOC1] 1
y[2,BOC2] 1
....Problem is that in every station there can be only installed one station. So one combination has to be 1 and all other 0, e.g. :
y[1,BTC1] 0
y[1,BTC2] 1
y[1,BOC1] 0
y[1,BOC2] 0To fix this issue that all are 1 I thought of a constraint like this:
m.addConstrs((quicksum(y[i,j] for i in stations) == 1 for j in types))
This wont work for me as it makes to optimization infeasible. How could I define a constraint that makes something like that? Thanks
0 -
The constraint makes sense for what you want to achieve. To find the issue of infeasibility, you should compute an IIS as described in How do I determine why my model is infeasible?
0 -
The constraints for the stations are doing what they should if I look at the LP file. What happens is, if I add the constraints to the model, the decision variable x becomes binary and generates no longer float numbers. Therefore the tasks dont fit in the stations because the capacity is exceeded and this is the reason why the model is infeasible.
It is very strange because why does this constraints effect the fractionality of x in any way? x should split up and take float numbers between 0 and 1 so that the Capacity of stations can be used fully. It should look like this for example:
x[BTC1,w01] 0
x[BTC1,w02] 5.6521739130434778e-01
x[BTC1,w03] 1
x[BOC1,w01] 1
x[BOC1,w02] 4.3478260869565222e-01
x[BOC1,w03] 0But when I add the constraints of:
m.addConstrs((quicksum(y[i,j] for i in stations) == 1 for j in types))
It rather looks like this:
x[BTC1,w01] 0
x[BTC1,w02] 0
x[BTC1,w03] 0
x[BOC1,w01] 1
x[BOC1,w02] 1
x[BOC1,w03] 1x only takes binary values which results in tasks not fitting in the number of stations and therefore it is infeasible. Remembering that the constraints for the stations are totally fine doing what they should
0 -
Did you try computing the IIS as proposed in my previous message?
Could you please post the current state of your model? I think it changed a lot with all the fixes you mentioned.
0 -
I did minimize the code for easier understanding and it looks like this:
stations = [1,2,3]
C = 200
C_min = 50
tasks = ['w01','w02','w03']
durations = {('BTC1','w01'):150, ('BTC1','w02'):115, ('BTC1','w03'):135,
('BOC1','w01'):150, ('BOC1','w02'):115, ('BOC1','w03'):135}
successors = {'w01':[],'w02':['w01'],'w03':['w02']}types = ['BTC1','BOC1']
f_cost = {(1,'BTC1'): 50000,(1,'BOC1'): 40000,
(2,'BTC1'): 50000,(2,'BOC1'): 40000,
(3,'BTC1'): 50000,(3,'BOC1'): 40000}
edges, v_cost = gp.multidict({
('BTC1','w01'): [15],
('BTC1','w02'): [30],
('BTC1','w03'): [45],
('BOC1','w01'): [50],
('BOC1','w02'): [80],
('BOC1','w03'): [110]
})x_arcs = [(i,j) for i in types for j in tasks]
y_arcs = [(i,j) for i in stations for j in types]m = gp.Model()
y = m.addVars(y_arcs, vtype = GRB.BINARY, name = 'y')
x = m.addVars(x_arcs, vtype = GRB.CONTINUOUS, name = 'x',ub=1)m.addConstrs(quicksum(x[i,j] for i in types) == 1 for j in tasks)
m.addConstrs(quicksum(durations[k,j]*x[k,j] for j in tasks) <= C*y[i,k] for i in stations for k in types)
#m.addConstrs(quicksum(y[i,j] for j in types) == 1 for i in stations)
m.addConstrs(quicksum(x[j,k] for j in types) <=
quicksum(x[j,i] for j in types)
for i in tasks for k in successors[i] if i!= 'w01')m.setObjective(x.prod(v_cost) + y.prod(f_cost), GRB.MINIMIZE)
m.optimize()The constraint markes as comment is the last one I added which constraints that every station has only 1 type of machine. Without it the code works fine. With it as I said the x variable is not Continuous between 0 and 1 anymore. When computing the IIS I get this:
Computing Irreducible Inconsistent Subsystem (IIS)... Constraints | Bounds | Runtime Min Max Guess | Min Max Guess | -------------------------------------------------------------------------- 0 14 - 0 12 - 0s 6 6 - 0 0 - 0s IIS computed: 6 constraints, 0 bounds IIS runtime: 0.02 seconds (0.00 work units)In the model.ilp File I get this:
Minimize
Subject To
R0: x[BTC1,w01] + x[BOC1,w01] = 1
R1: x[BTC1,w02] + x[BOC1,w02] = 1
R2: x[BTC1,w03] + x[BOC1,w03] = 1
R3: - 200 y[1,BTC1] + 150 x[BTC1,w01] + 115 x[BTC1,w02] + 135 x[BTC1,w03]
<= 0
R4: - 200 y[1,BOC1] + 150 x[BOC1,w01] + 115 x[BOC1,w02] + 135 x[BOC1,w03]
<= 0
R9: y[1,BTC1] + y[1,BOC1] = 1
Bounds
x[BTC1,w01] free
x[BTC1,w02] free
x[BTC1,w03] free
x[BOC1,w01] free
x[BOC1,w02] free
x[BOC1,w03] free
Binaries
y[1,BTC1] y[1,BOC1]
End0 -
Note that the default lower bounds value for variables in Gurobi is \(0\), so your \(x\) variables are all between \(0\leq x \leq 1\). From the IIS, you can deduce what is going wrong. Let us set \(y_{1,BTC1} = 1\). Then \(y_{1,BOC1}=0\). With \(y_{1,BOC1}=0\), from
\[\begin{align*}
150 x_{BOC1,w01} + 115 x_{BOC1,w02} + 135 x_{BOC1,w03} \leq 200 \cdot 0
\end{align*}\]follows \(x_{BOC1,w01} = x_{BOC1,w02} = x_{BOC1,w03} = 0\). From constraints \(R0,R1,R2\), we get
\(x_{BTC1,w01} = x_{BTC1,w02} = x_{BTC1,w03} = 1\) which leads to constraint \(R3\)\[\begin{align*}
150\cdot 1+ 115 \cdot 1 + 135 \cdot 1 \leq 200 \cdot 1
\end{align*}\]which is infeasible. An analogous statement happens for \(y_{1,BOC1}=1\).
One possibility to make your model feasible would be to allow for negative \(x\) values. However, I don't know the details of your model and whether this would make any sense. Alternatively, you could increase the coefficient \(200\) to something bigger or adjust the constraints \(R0,R1,R2\).
Best regards,
Jaromił0 -
Hi Jaromił Najman I just wanted to inform you that all the problems I had are solved. The program code is doing all things I want to and the constraints are also on point, and at the same time the solutions are logic and sensible. I wanted to say huuuuge thanks to you, you supported me within these few days with my, as for my levels, very complex issue. I lastly want to plot the results, if I get any problems with that I will open a new Question and define this one as closed. I wish you a nice weekend and god bless you!
0 -
Hi Harun Gül,
Was just going through your problem. I have similar issues as yours. I was wondering if we can connect to discuss the precedence of task code and how did you manage to resolve it,Best,
Roy
0
Please sign in to leave a comment.
Comments
19 comments