5

I will try to keep my question short and simple. If you need any further information, please let me know.

I have an MIP, implemented in Python with the package PuLP. (Roughly 100 variables and constraints) The mathematical formulation of the problem is from a research paper. This paper also includes a numerical study. However, my results differ from the results of the authors.

My problem variable is called prob

prob = LpProblem("Replenishment_Policy", LpMinimize)

I solve the problem with prob.solve() LpStatus returns Optimal

When I add some of the optimal (paper) results as contraints, I get a slightly better objective value. Same goes for constraining the objecive function to a slightly lower value. The LpStatus remains Optimal.

original objective value: total = 1704.20
decision variable: stock[1] = 370

adding constraints: prob += stock[1] == 379
new objective value: 1704.09

adding constraints: prob += prob.objective <= 1704
new objective value: 1702.81

My assumption is that PuLP's solver approximates the solution. The calculation is very fast, but apparently not very accurate. Is there a way I can improve the accuracy of the solver PuLP is using? I am looking for something like: prob.solve(accuracy=100%). I had a look at the documentation but couldn't figure out what to do. Are there any thoughts what the problem could be?

Any help is appreciated. Thanks.

Axel
  • 2,545
  • 2
  • 18
  • 30
  • 4
    You are probably using the default solver, CBC. In that case, I use `prob.solve(solvers.PULP_CBC_CMD(fracGap=0.01))` where fracGap is the tolerance for the optimal solution (within 1% accuracy). If you decrease this your accuracy will increase. I don't think passing 0 would work because of the floating point issues but of course you can decrease it a lot for an LP. – ayhan Jun 28 '17 at 17:59
  • 3
    ```The calculation is very fast, but apparently not very accurate``` This indicates a problem or a misunderstanding, as the difference between a valid solution using the default config and the optimal solution should be less than 0.1%. – sascha Jun 28 '17 at 18:03
  • Thanks for the quick reply. Unfortunately @ayhan's suggestion did not improve the results. As @sascha mentioned correctly, the difference between the two valid solutions is larger than a default fracGap. This probably means the error lies somewhere else. Does the solver's `fracGap` also apply for constraints? If not, how/where could I specify it? – Axel Jun 28 '17 at 19:22
  • The fracGap applies to solutions. All solutions have to be feasible = constraints are fulfilled! Make sure your solution is a valid solution, meaning: check the status of the solver. – sascha Jun 28 '17 at 19:31
  • @Axel I doubt that it is the solver's accuracy that is the cause of the difference. Can you post the obj function value that you are getting vs. the published one? I suspect one of the coefficients may be slightly off. – Ram Narasimhan Jun 29 '17 at 04:21
  • Unfortunately the paper, does not state the optimal value, but only the results of the variables. In the meanwhile, I have found the issue of my problem and posted an answer which will hopefully elaborate and help future visitors at the same time. – Axel Jun 29 '17 at 16:25

1 Answers1

2

The answer to my question was given by ayhan: To specify the accuracy of the solver, you can use the fracGap argument of the selected solver.

prob.solve(solvers.PULP_CBC_CMD(fracGap=0.01))

However, the question I asked, was not aligned with the problem I had. The deviation of the results was indeed not a matter of accuracy of the solver (as sascha pointed out in the comments).

The cause to my problem: The algorithm I implemented was the optimization of the order policy parameters for a (Rn, Sn) policy under non-stationary, stochastic demand. The above mentioned paper is: Tarim, S. A., & Kingsman, B. G. (2006). Modelling and computing (R n, S n) policies for inventory systems with non-stationary stochastic demand. European Journal of Operational Research, 174(1), 581-599.

The algorithm has two binary variables delta[t] and P[t][j]. The following two constraints only allow values of 0 and 1 for P[t][j], as long as delta[t] is defined as a binary.

for t in range(1, T+1):
    prob += sum([P[t][j] for j in range(1, t+1)]) == 1

    for j in range(1, t+1):
        prob += P[t][j] >= delta[t-j+1] - sum([delta[k] for k in range(t-j+2, t+1)])

Since P[t][j] can only take values of 0 or 1, hence being a binary variable, I declared it as follows:

for t in range(1, T+1):
    for j in range(1, T+1):
        P[t][j] = LpVariable(name="P_"+str(t)+"_"+str(j), lowBound=0, upBound=1, cat="Integer")

The objective value for the minimization returns: 1704.20

After researching for a solution for quite a while, I noticed a part of the paper that says:

... it follows that P_tj must still take a binary value even if it is declared as a continuous variable. Therefore, the total number of binary variables reduces to the total number of periods, N.

Therefore I changed the cat argument of the P[t][j] variable to cat="Continuous". Whithout changing anything else, I got the lower objective value of 1702.81. The status of the result shows in both cases: Optimal

I am still not sure how all these aspects are interrelated, but I guess for me this tweek worked. For everyone else who is directed to this question, will probably find the necessary help with the answer given at the top of this post.

Axel
  • 2,545
  • 2
  • 18
  • 30