The issue to writing the multiple solutions

[abb-omidi]

Dear support team,

I am working on the column generation procedure to solve a scheduling problem. For that, as the sub-problem is a knapsack form, I would like to count all of the extreme points from the pricing and add those entirely to the master problem. I saw the related questions, but I am unsure whether the issue of taking the multiple optimal solutions is being fixed. What I am trying is:

mdl = Model()

# Defining the variables
# Defining the constraints

mdl.setPresolve(SCIP_PARAMSETTING.OFF)
mdl.setHeuristics(SCIP_PARAMSETTING.OFF)
mdl.setSeparating(SCIP_PARAMSETTING.OFF)
mdl.setBoolParam("constraints/countsols/collect", True)
mdl.setLongintParam("constraints/countsols/sollimit", 20)

mdl.writeProblem("Test.lp")

mdl.count()
nsols = mdl.getNCountedSols()
sols = mdl.getSols()

The result is an empty list!!!

By solving the problem with SCIP, I can get the results as following:

# | y_2 | y_1 | x_item2_2 | x_item2_1 | x_item1_2 | x_item1_1
-- | -- | -- | -- | -- | -- | --
1(1) | 0 | 0 | 0 | 0 | 0 | 0
2(2) | 0 | 1 | 0 | 0 | 0 | 0
2(3) | 0 | 1 | 0 | 1 | 0 | 0
3(4) | 0 | 1 | 0 | 0 | 0 | 1
4(5) | 1 | 1 | 0 | 0 | 0 | 1
4(6) | 1 | 1 | 1 | 0 | 0 | 1
5(7) | 1 | 1 | 0 | 0 | 1 | 1
6(8) | 1 | 0 | 0 | 0 | 0 | 0
6(9) | 1 | 0 | 1 | 0 | 0 | 0
7(10) | 1 | 1 | 0 | 0 | 0 | 0
7(11) | 1 | 1 | 0 | 1 | 0 | 0
7(12) | 1 | 1 | 1 | 0 | 0 | 0
7(13) | 1 | 1 | 1 | 1 | 0 | 0
8(14) | 1 | 0 | 0 | 0 | 1 | 0
9(15) | 1 | 1 | 0 | 0 | 1 | 0
9(16) | 1 | 1 | 0 | 1 | 1 | 0

However, I am trying to use the method mdl.writeSol(), but I am not aware of how to define the required argument solution to test and see what's really happened. I was wondering if, you could take a look and hint to me about that.

Installing collected packages: pyscipopt
Successfully installed pyscipopt-5.1.1
Install by PIP

All the best
Abbas

[Joao-Dionisio]

Hey @abb-omidi! I think the issue is with the count method. Look at issue #838, it might be helpful to you.

In regards to column generation, I think using something like count will not be very helpful. You only want to get the negative reduced cost columns, and count will try to find every solution. Sure you can add a constraint to exclude the non-negative reduced cost ones, but even still, usually just keeping the ones found by SCIP when optimizing is enough. You want to get the updated duals for the next iteration, spending extra time finding suboptimal solutions does not seem the best approach.

[abb-omidi]

Dear @Joao-Dionisio,

Thanks for your comments. Also, the new online library is very useful.
I think the count method works well as it can truly gather feasible solutions. What I doubt works fine comes back to the writing function write allsolutions, which in SCIP writes the solution found in an external file. While it seems there is no such a method in the PySCIPOpt. The only defined method is writeSol which can display one of the gathering solutions in the best case. B.T.W I cannot find a way to write all of the solutions. However, having a nice function like count would be useful in many cases.

I tried to write an event handler to see what happened for BESTSOLFOUND as long as I called count instead of optimize, but it threw an error on the execution function:

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
[<ipython-input-8-1f83dbe54bc3>](https://localhost:8080/#) in <cell line: 77>()
     76 
     77 if __name__ == "__main__":
---> 78   run_event()

[<ipython-input-8-1f83dbe54bc3>](https://localhost:8080/#) in run_event()
     70   print(calls)
     71   assert 'eventinit' in calls
---> 72   assert 'eventexec' in calls
     73    assert 'eventexit' in calls
     74    assert len(calls) == 3

AssertionError:

Regarding the column generation, in the first step, I have tried to generate the polyhedron from scratch to see the behavior of the master problem objective bounds Vs the objective of the compact formulation when the columns are added entirely at the beginning. In the second step, I would go with what you proposed. (That I will get back to you soon)

All the best

[Joao-Dionisio]

This is what @Opt-Mucca (who was also the one who made the new documentation, very useful indeed!) had to say about it:

If you enable countsols then SCIP is essentially hiding the solutions from itself. So it stores the solutions internally and gives you a valid count of solutions found, but SCIP itself does not have access to them (to do things like prune the search tree, etc, or terminate the search early).

I agree it would be nice if count worked the same way optimize does, I just don't know how big a change this would end up being.

--

As for your error in the callback, I think it might make sense, but it's difficult to be sure without the full code. From what I could understand of count, it basically throws away the objective function and it simply becomes a feasibility problem. Every time it finds a new solution, it adds a cut removing it, and reoptimizes. Without the objective function, maybe SCIP doesn't update BESTSOLFOUND. You're saying that the same callback works when using optimize instead of count?

[abb-omidi]

It would be somewhat strange if SCIP could write the feasible solutions itself without throwing any issues, but PySCIPOpt cannot do that! It would be great if someone could take a closer look at the mapping functions between SCIP and PySCIPOpt to fix that.

An MRE to reproduce the error and the issue would be:

from pyscipopt import Model, Eventhdlr, SCIP_RESULT, SCIP_EVENTTYPE, SCIP_PARAMSETTING, quicksum

class MyEvent(Eventhdlr):

    def __init__(self):
        pass

    def eventinit(self):
        self.model.catchEvent(SCIP_EVENTTYPE.BESTSOLFOUND, self)

    def eventexit(self):
        self.model.dropEvent(SCIP_EVENTTYPE.BESTSOLFOUND, self)

    def eventexec(self, event):
        assert event.getType() == SCIP_EVENTTYPE.BESTSOLFOUND
        print("best solution found")
        result = SCIP_RESULT.FOUNDSOL
    
    
def model_opt():
    s = Model()
    s.redirectOutput()
    s.printVersion()
    s.setBoolParam("constraints/countsols/collect", True)
    s.setLongintParam("constraints/countsols/sollimit", 20)
    s.setPresolve(SCIP_PARAMSETTING.OFF)
    s.setSeparating(SCIP_PARAMSETTING.OFF)
    s.setHeuristics(SCIP_PARAMSETTING.OFF)
    print("--------------------------------------")

    eventhdlr = MyEvent()
    s.includeEventhdlr(eventhdlr, "BESTSOLFOUND", "python event handler to catch BESTSOLFOUND")
    eventhdlr.model = s
    return s

def run_event():
    s = model_opt()

    x_item1_1 = s.addVar("x_item1_1", vtype="B", obj=0)
    x_item2_1 = s.addVar("x_item2_1", vtype="B", obj=0)
    x_item1_2 = s.addVar("x_item1_2", vtype="B", obj=0)
    x_item2_2 = s.addVar("x_item2_2", vtype="B", obj=0)
    y_1 = s.addVar("y_1", vtype="B", obj=0)
    y_2 = s.addVar("y_2", vtype="B", obj=0)

    s.addCons(+3 * x_item1_1 +5 * x_item2_1 -5 * y_1 <= +0)
    s.addCons(+3 * x_item1_2 +5 * x_item2_2 -5 * y_2 <= +0)

    s.count()

    print()
    nsols = s.getNCountedSols()
    print("\nNumber of solution found: ", nsols)
    s.writeSol(s.count(), filename="origprob.sol", write_zeros=True)

    assert 'eventinit' in calls
    # assert 'eventexit' in calls
    assert 'eventexec' in calls

if __name__ == "__main__":
    run_event()

The result should be:

SCIP version 9.1.0 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver: Soplex 7.1.0] [GitHash: NoGitInfo]
Copyright (c) 2002-2024 Zuse Institute Berlin (ZIB)
--------------------------------------
counting forces parameter <misc/usesymmetry> to 0.
presolving:
presolving (0 rounds: 0 fast, 0 medium, 0 exhaustive):
 0 deleted vars, 0 deleted constraints, 0 added constraints, 0 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
 0 implications, 0 cliques
presolved problem has 6 variables (6 bin, 0 int, 0 impl, 0 cont) and 2 constraints
      2 constraints of type <linear>
transformed objective value is always integral (scale: 1)
Presolving Time: 0.00

a non-trivial solution comes in over <SCIP_DECL_CONSCHECK(consCheckCountsols)>; currently these solutions are ignored.
 time | node  | left  |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr|  dualbound   | primalbound  |  gap   | compl. 
  0.0s|     1 |     0 |     0 |     - |   601k |   0 |   6 |   2 |   2 |   0 |  0 |   0 |   0 | 0.000000e+00 |      --      |    Inf | unknown
  0.0s|     1 |     2 |     0 |     - |   601k |   0 |   6 |   2 |   2 |   0 |  1 |   0 |   0 | 0.000000e+00 |      --      |    Inf | unknown

SCIP Status        : problem is solved [infeasible]
Solving Time (sec) : 0.01
Solving Nodes      : 21
Primal Bound       : +1.00000000000000e+20 (0 solutions)
Dual Bound         : +1.00000000000000e+20
Gap                : 0.00 %


Number of solution found:  16
objective value:                                    0
x_item1_1                                           0 	(obj:0)
x_item2_1                                          -0 	(obj:0)
x_item1_2                                           1 	(obj:0)
x_item2_2                                           0 	(obj:0)
y_1                                                 1 	(obj:0)
y_2                                                 1 	(obj:0)
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
[<ipython-input-20-2bcbd478777f>](https://localhost:8080/#) in <cell line: 61>()
     60 
     61 if __name__ == "__main__":
---> 62     run_event()

[<ipython-input-20-2bcbd478777f>](https://localhost:8080/#) in run_event()
     57     assert 'eventinit' in calls
     58     # assert 'eventexit' in calls
---> 59     assert 'eventexec' in calls
     60 
     61 if __name__ == "__main__":

AssertionError:

I was wondering if, you could take a look at this and see if it can be extended to write the rest of the solution.

Regards

[abb-omidi]

Just as a follow-up question:

would you say please, how can we relax an integer problem without changing the variables type?
As I could find in the examples, there are some functions like relaxed or LP, but I am unsure if I understand how to use that in my case.

Also. in the reference you provided for B&P on the bin-packing problem, it seems it can be done by getTransformedCons, getDualfarkasLinear, and getDualsolLinear. I am wondering if, is there any simple example to use these methods to fully understand how they work.

P.S: @Joao-Dionisio, I could resolve some of my issues regarding relaxation and duals. Pls, let me check something, and I will get back to you asap.

All the best

[abb-omidi]

Dear support team,

Is it a necessery thing to define a pricer to capture the variables attributes to calculate the appropriate added columns or we can still use a loop to iterate over the extreme points/rays and update pricing problem and add column inside that?

Best

[Joao-Dionisio]

how can we relax an integer problem without changing the variables type?

there is the relax() method.

Also. in the reference you provided for B&P on the bin-packing problem, it seems it can be done by getTransformedCons, getDualfarkasLinear, and getDualsolLinear. I am wondering if, is there any simple example to use these methods to fully understand how they work.

There is not, you can take a look at the documentation of the functions, which should explain them well enough.

Is it a necessery thing to define a pricer to capture the variables attributes to calculate the appropriate added columns or we can still use a loop to iterate over the extreme points/rays and update pricing problem and add column inside that?

There is nothing stopping you from using SCIP as an LP solver and do the column-generation process yourself, it's just extra work.

[Opt-Mucca]

We will work on adding a function for getting the counted solutions. The issue is that they are stored in the constraint handler, not SCIPs traditional solution storage. We just need to write an appropriate function.

[abb-omidi]

Dear @Opt-Mucca,
Thank you so much for your interest. Please, inform us when it is down. (🙏)

Dear @Joao-Dionisio,
Thanks for your answers and comments on how to run the code on GitHub. Also, I have tried to write the code I proposed in my last email in a modeling language. It works fine, but it cannot guarantee that it will reach optimality as it does not contain any tree search on the fly like the branching rule. (I actually implement price&branch). However, I need much time to adapt it to PySCIPOpt, specifically from the data structure point of view.

I am currently working on the BnP file you provided and I have some sort of questions. In order to implement branch&price and as far as I know,

1- we need to formulate an extended master problem.
2- defining the pricing model.
3- Initiating the master to ignite the process.
4- Looping over the master and pricing.
5- Checking if the master solution is integral, otherwise, branching.

Now, my questions are:

• Would you confirm that my above intuition is correct?
• Is there any reason to separate the pricer from the knapsack model?
• I am not well familiar with Ryan-Foster, but would you confirm that adding the extra constraints in the pricer causes the knapsack polytope being unbounded? (I ask because in the code we obtain Farkas certificates!)
• Does the function solve_knapsack() solve towards the optimization process? If so, what is the role of the function solve_knapsack_with_constraints? (AFAIK, the latter should be used?)
• Is there a way to write the pricing problem as well as the master problem in each iteration? (I really would like to investigate the process step by step by hand and compare both results) Also, I think I have to set the node limit parameter to run iteratively.

All the best

[Joao-Dionisio]

Hey @abb-omidi .

1 - yes, your general understanding is correct.
2 - just to make the code tidier.
3 - I do not confirm that, a polytope cannot become unbounded by adding constraints. The Farkas certificates are mostly for the infeasibilities that might occur due to branching. They can also be useful for initializing column generation (see Farkas pricing)
4 - the knapsack without constraints is useful only for column generation at the root node, you are correct that the latter should be used.
5 - yes, you can use writeProblem() at the appropriate places.

[abb-omidi]

Dear @Joao-Dionisio,

Would you confirm that:

• The written problem of BnP in Github and the provided one before taking place co@work work the same? It seems the Github version in some classes has an entirely different data structure rather than the previous one, however as I checked in some instances the provided results of both are the same.
• What actually the Farkas pricing works? (As I cannot find any appropriate documents.)
• In the x[i] = m.addVar(vtype="B", name=f"x{i}", obj=values[i]) at knapsack class, how values[i] may take the dual value of the master constraints?
• I have modified the code so that it can write the problems in each iteration of the algorithm. It seems that the master problem is only written at the first and the end of runtime. Is this normal behavior? I expected the master to be called many times in the solving process to update the duals and the rest of the needed parameters.
• What is the difference between result[0] in the min_red_cost = 1 - result[0] and result[1] that be returned by function pricing_solver?
• Is it possible to turn the heuristics off as long as the master and pricing problems are being solved?

All the best

[abb-omidi]

Dear @Joao-Dionisio,

I think there is a strange behavior in the running process of the BnP code. (Still In both versions as I said in my previous comment). The issue is that as long as the size of the items has exceeded bins capacity it can produce an infeasible solution as a feasible one. E.g. by the following simple data:

capacity = 5
sizes = [7,7,7]

the solver shows the following output while the compact model is running without any issues.

SCIP version 9.1.0 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver: Soplex 7.1.0] [GitHash: NoGitInfo]
Copyright (c) 2002-2024 Zuse Institute Berlin (ZIB)
presolving:
presolving (0 rounds: 0 fast, 0 medium, 0 exhaustive):
 0 deleted vars, 0 deleted constraints, 0 added constraints, 0 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
 0 implications, 0 cliques
presolved problem has 3 variables (3 bin, 0 int, 0 impl, 0 cont) and 3 constraints
      3 constraints of type <linear>
Presolving Time: 0.00

 time | node  | left  |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr|  dualbound   | primalbound  |  gap   | compl. 
* 0.0s|     1 |     0 |     3 |     - |    LP  |   0 |   3 |   3 |   3 |   0 |  0 |   0 |   0 |      --      | 3.000000e+00 |    Inf | unknown
  0.0s|     1 |     0 |     3 |     - |   591k |   0 |   3 |   3 |   3 |   0 |  0 |   0 |   0 | 3.000000e+00 | 3.000000e+00 |   0.00%| unknown
  0.0s|     1 |     0 |     3 |     - |   591k |   0 |   3 |   3 |   3 |   0 |  0 |   0 |   0 | 3.000000e+00 | 3.000000e+00 |   0.00%| unknown

SCIP Status        : problem is solved [optimal solution found]
Solving Time (sec) : 0.02
Solving Nodes      : 1
Primal Bound       : +3.00000000000000e+00 (1 solutions)
Dual Bound         : +3.00000000000000e+00
Gap                : 0.00 %

objective value 3.0
Solution:
x_[0], x_[1], x_[2], 

solving compact model: 
SCIP version 9.1.0 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver: Soplex 7.1.0] [GitHash: NoGitInfo]
Copyright (c) 2002-2024 Zuse Institute Berlin (ZIB)
presolving:
presolving (1 rounds: 1 fast, 0 medium, 0 exhaustive):
 9 deleted vars, 3 deleted constraints, 0 added constraints, 9 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
 0 implications, 3 cliques
presolving detected infeasibility
Presolving Time: 0.00

SCIP Status        : problem is solved [infeasible]
Solving Time (sec) : 0.02
Solving Nodes      : 0
Primal Bound       : +1.00000000000000e+20 (0 solutions)
Dual Bound         : +1.00000000000000e+20
Gap                : 0.00 %

I tried to write the solve_knapsack_with_constraints problem and the result shows a somewhat irregular model.

\ SCIP STATISTICS
\   Problem name     : t_knapsack_with_constraints
\   Variables        : 3 (3 binary, 0 integer, 0 implicit integer, 0 continuous)
\   Constraints      : 0
Maximize
 Obj: +1 t_x0 +1 t_x1 +1 t_x2
Subject to
Bounds
 0 <= t_x0 <= 0
 0 <= t_x1 <= 0
 0 <= t_x2 <= 0
Binaries
 t_x0 t_x1 t_x2
End

Would you say please, it is a bug or I am missing something?
(I think we need to add an event handler to check whether the subproblem is feasible or not. If not, return an infeasible flag and stop the process.).

Best regards

[Joao-Dionisio]

• They should work more or less the same, but the co@work repo is better tested, so you should stick to that one.
• Section 4.1.2 of Marco Lubbecke's PhD thesis "Generic Branch-Cut-and-Price"
• the price of the items in the knapsack subproblem is given by the duals of the master, you can take a look at marco's tutorial online
• yes, this is normal behavior. SCIP is working more efficiently than just creating the master again and again at each iteration
• 1 - result[0] is the reduced cost. reduced[1] is the pattern, the actual solution, if I recall correctly.
• You mean turning them off at the beginning? Yes, it is possible with model.setHeuristics(0), if my memory doesn't fail me.

That seems to be a bug, Abbas, thank you. Please don't open any more issues about co@work in this repo, do so in the co@work repo, please. I'll look later into what you talked about.

[abb-omidi]

Dear @Joao-Dionisio,

Many thanks for your answers, specifically for the 4th bullet.
Just for clarification would you say you mean by reduced[1] is actually result[1]?

Best regards