Adding primal solutions of original problem after the model has been transformed

[SomovMike]

Hello! I try to add my own heuristic, which called while branching and give primal solutions to SCIP. This heuristic works only with original problem, but after presolving SCIP add new different variables and constraints. I've tried to create solution with CreateSol() method, and then put solution to SCIP with trySol() or addSol(), but as I said my heuristic get only values of original variables, and CreateSol() creates solution with transformed variables. I guess trySol() and addSol() also require transformed variables. Is it possible to put original solution while solving? Only way I see to off presolving, but it is not the best option for me.
Thank you in advance!

[Joao-Dionisio]

Hello, @SomovMike ! I haven't used SCIP during branching much, but I can take a look. Can you provide a code example of the problem you're having?

I don't know if it works, but there is also the method createPartialSol(), where you only provide the values to a subset of variables.

EDIT: Maybe you have already seen these, but here is a link to adding Primal Heuristics in SCIP, and a code example of adding a heuristic in PySCIPOpt.

EDIT2: Does this question in StackOverflow help? Link

[SomovMike]

Thank you for your quick answer! But I do not understand how createPartialSol() should work, this function also create solution with all transformed variables, when I put only values of original variables in it and try to add this solution with addSol() function it always return False, as I understood it means that this solution didn't get new primal bound

[SomovMike]

Here is example of output of createPartialSol() function after presolving, I can provide only values of "t_x_i" variables, but as you can see there are a lot of new binreform variables.
{'t_x_0': 1e+98, 't_x_1': 1e+98, 't_x_2': 1e+98, 't_x_3': 1e+98, 't_x_4': 1e+98, 't_x_5': 1e+98, 't_x_6': 1e+98, 't_x_7': 1e+98, 't_x_8': 1e+98, 't_x_9': 1e+98, 't_x_10': 1e+98, 't_x_11': 1e+98, 't_x_12': 1e+98, 't_x_13': 1e+98, 't_x_14': 1e+98, 't_x_15': 1e+98, 't_x_16': 1e+98, 't_x_17': 1e+98, 't_x_18': 1e+98, 't_x_19': 1e+98, 't_x_20': 1e+98, 't_x_21': 1e+98, 't_x_22': 1e+98, 't_x_23': 1e+98, 't_x_24': 1e+98, 't_x_25': 1e+98, 't_x_26': 1e+98, 't_x_27': 1e+98, 't_x_28': 1e+98, 't_x_29': 1e+98, 't_x_30': 1e+98, 't_x_31': 1e+98, 't_x_32': 1e+98, 't_x_33': 1e+98, 't_x_34': 1e+98, 't_x_35': 1e+98, 't_x_36': 1e+98, 't_x_37': 1e+98, 't_x_38': 1e+98, 't_x_39': 1e+98, 't_x_40': 1e+98, 't_x_41': 1e+98, 't_x_42': 1e+98, 't_x_43': 1e+98, 't_x_44': 1e+98, 't_x_45': 1e+98, 't_x_46': 1e+98, 't_x_47': 1e+98, 't_x_48': 1e+98, 't_x_49': 1e+98, 't_x_50': 1e+98, 't_x_51': 1e+98, 't_x_52': 1e+98, 't_x_53': 1e+98, 't_x_54': 1e+98, 't_x_55': 1e+98, 't_x_56': 1e+98, 't_x_57': 1e+98, 't_x_58': 1e+98, 't_x_59': 1e+98, 't_x_60': 1e+98, 't_x_61': 1e+98, 't_x_62': 1e+98, 't_x_63': 1e+98, 'binreform_t_x_1_t_x_8': 0.0, 'binreform_t_x_7_t_x_8': 0.0, 'binreform_t_x_0_t_x_9': 0.0, 'binreform_t_x_2_t_x_9': 0.0, 'binreform_t_x_1_t_x_10': 0.0, 'binreform_t_x_3_t_x_10': 0.0, 'binreform_t_x_2_t_x_11': 0.0, 'binreform_t_x_4_t_x_11': 0.0, 'binreform_t_x_3_t_x_12': 0.0, 'binreform_t_x_5_t_x_12': 0.0, 'binreform_t_x_4_t_x_13': 0.0, 'binreform_t_x_6_t_x_13': 0.0, 'binreform_t_x_5_t_x_14': 0.0, 'binreform_t_x_7_t_x_14': 0.0, 'binreform_t_x_0_t_x_15': 0.0, 'binreform_t_x_6_t_x_15': 0.0, 'binreform_t_x_1_t_x_16': 0.0, 'binreform_t_x_7_t_x_16': 0.0, 'binreform_t_x_9_t_x_16': 0.0, 'binreform_t_x_15_t_x_16': 0.0, 'binreform_t_x_0_t_x_17': 0.0, 'binreform_t_x_2_t_x_17': 0.0, 'binreform_t_x_8_t_x_17': 0.0, 'binreform_t_x_10_t_x_17': 0.0, 'binreform_t_x_1_t_x_18': 0.0, 'binreform_t_x_3_t_x_18': 0.0, .....}

[Joao-Dionisio]

My suggestion with createPartialSol() was that maybe by assigning values to the "t_x_i" variables, the remaining variables would be easily found by SCIP (I have no idea if that is the case, though!). It's strange that you are having problems with createPartialSol() this way. Can you provide a small example of the code you're running, so I can try to reproduce it, please?

In the link to StackOverflow I sent above, Gregor says the following:

You can always query solution values of original variables, even from transformed solutions, using SCIPgetSolVal()

So it seems that SCIP somehow saves the original problem. Looking at the documentation, I found SCIPgetOrigVars(), that gets the original variables, but I don't think that this method has been brought over to PySCIPOpt yet.

[SomovMike]

from pyscipopt import Model, quicksum, Heur, SCIP_HEURTIMING, SCIP_RESULT import numpy as np

Create Heuristic which get optimal solution of original variables

class MyHeur(Heur):
    
    def __init__(self, optimal_solution):
        self.optimal_solution = optimal_solution
    
    
    def heurexec(self, heurtiming, nodeinfeasible):
        scip_sol = self.model.createPartialSol(self)
        for var in self.model.getVars()[:-1]:
            scip_sol[var] = self.optimal_solution[var.getIndex()]
        accepted = self.model.addSol(scip_sol)
        
        if accepted:
            print("Accepted")
            return {"result": SCIP_RESULT.FOUNDSOL}
        else:
            print("Not accepted")
            return {"result": SCIP_RESULT.DIDNOTFIND}

Function to create TSP model

def create_TSP_model(c, V):
    model = Model("tsp")

    x = {}
    for i in V:
        for p in V:
            x[i,p] = model.addVar(vtype="B", name="x(%s,%s)"%(i,p))

    obj_var = model.addVar(vtype="C", name="obj_var")

    for i in V:
        model.addCons(quicksum(x[i, p] for p in V) == 1)
    for p in V:
        model.addCons(quicksum(x[i, p] for i in V) == 1)

    tmp = [i for i in V] + [0]
    model.addCons(obj_var >= quicksum(c[i, j]*quicksum(x[i, tmp[p]]*x[j, tmp[p+1]] for p in V) for i in V for j in V))

    model.setObjective(obj_var, "minimize")
    
    return model

Data

c = np.array([[ 0.  , 64.33, 79.38, 41.15, 42.2 , 78.6 , 49.48, 46.17],
       [64.33,  0.  , 40.8 , 40.46, 47.8 , 77.67, 79.51, 33.02],
       [79.38, 40.8 ,  0.  , 38.83, 41.01, 46.4 , 66.73, 33.24],
       [41.15, 40.46, 38.83,  0.  ,  7.62, 46.1 , 39.05,  8.06],
       [42.2 , 47.8 , 41.01,  7.62,  0.  , 40.16, 31.83, 14.87],
       [78.6 , 77.67, 46.4 , 46.1 , 40.16,  0.  , 38.08, 49.09],
       [49.48, 79.51, 66.73, 39.05, 31.83, 38.08,  0.  , 46.69],
       [46.17, 33.02, 33.24,  8.06, 14.87, 49.09, 46.69,  0.  ]])
V = range(c.shape[0])

Solve model without my heuristic and get optimal soltuion

model_without_heur = create_TSP_model(c, V)
model_without_heur.optimize()
variables = []
for var in model_without_heur.getVars():
    if var.name != "obj_var":
        variables.append(var)
optimal_solution = np.array([model_without_heur.getVal(var) for var in variables])

Then solve model with my heuristic and put optimal solution in it

model_with_heur = create_TSP_model(c, V)
heuristic = MyHeur(optimal_solution)
model_with_heur.includeHeur(heuristic, "PyHeur", "custom heuristic", "Y",
                              timingmask=SCIP_HEURTIMING.BEFORENODE, freq=1, freqofs=0, maxdepth=2)

Try to optimize this model

model_with_heur.optimize()

So I always got "Not accepted" message, but solution is optimal, I guess createPartialSol() also requires values of other variables, or maybe "binreform..." variables are kind of special.

[Joao-Dionisio]

I've been looking at the code, and I don't think it's a problem with the transformed variables. I changed the code a bit, to reflect what I have encountered around SCIP.

I also made a change here, to create a dictionary that maps the variable names to their values in the optimal solution. I find it easier to work with.

model_without_heur = create_TSP_model(c, V)
model_without_heur.optimize()
optimal_solution = {}
for var in model_without_heur.getVars():
    optimal_solution[var.name] = model_without_heur.getVal(var)

model_with_heur = create_TSP_model(c, V)
heuristic = MyHeur(optimal_solution,model_without_heur) # I also pass the original model to the heuristic, just for testing

And the heuristic:

class MyHeur(Heur):
    
    def __init__(self, optimal_solution,model_without_heur):
        self.optimal_solution = optimal_solution
        self.model_without_heur = model_without_heur
    
    def heurexec(self, heurtiming, nodeinfeasible):
        scip_sol = self.model.createPartialSol()
        for var in self.model.getVars()[:-1]:
            self.model.setSolVal(scip_sol,var,self.optimal_solution[var.name]) # I think this is the standard way of assigning values to a solution 

        model_var_names = [var.name for var in self.model.getVars()] 
        model_without_heur_var_names = [var.name for var in self.model_without_heur.getVars()]

        print(model_var_names == model_without_heur_var_names) # This prints True, so they have the same variable names

        accepted = self.model.addSol(scip_sol)
        
        if accepted:
            print("Accepted")
            return {"result": SCIP_RESULT.FOUNDSOL}
        else:
            print("Not accepted")
            return {"result": SCIP_RESULT.DIDNOTFIND}

So, at this moment, the variables in the current model and the model without the heuristic have the same names, but I think you are right that accepted should be True as well. I will take a better look at this, but at the moment I still haven't found a solution.

EDIT: If you change the code for the heuristic that I wrote above so that you have scip_sol = self.model.createSol() instead of scip_sol = self.model.createPartialSol(), then accepted becomes True. However, the resulting solution is detected as infeasible, for some reason.

I think I might be getting out of my league here...

[SomovMike]

Thank you very much for for your help! Maybe I can somehow get values of transformed variables from original variables, and then create full solution with CreateSol() function, but I can't find any function of method that can do that:(

[Joao-Dionisio]

I think I can already do that, but for some reason SCIP says that the resulting solution is not feasible

It says that the solution is accepted. If I use trySol (which checks for feasibility), then he says that the solution is not feasible:

Maybe the problem is the way SCIP is trying to deal with your obj_var variable, due to the nonlinear objective function? The infeasibilities are indeed with the transformed variables, so maybe you're right... either way, I don't think I know how to solve your problem, sorry...

[SomovMike]

EDIT: If you change the code for the heuristic that I wrote above so that you have scip_sol = self.model.createSol() instead of scip_sol = self.model.createPartialSol(), then accepted becomes True. However, the resulting solution is detected as infeasible, for some reason.

As I understood solution infeasible because of additional variables, they just have default 0 value when you create solution with CreateSol(), but it is really strange that this solution is accepted. Feasible or infeasible, doesn't matter. It will be good for me even it just speed up solution process.

[SomovMike]

Maybe the problem is the way SCIP is trying to deal with your obj_var variable, due to the nonlinear objective function? The infeasibilities are indeed with the transformed variables, so maybe you're right... either way, I don't think I know how to solve your problem, sorry...

Thank you anyway!

[SomovMike]

Found interesting thing, if we rewrite model like this, just change "obj_var >= ..." on "obj_var + 1 >= ..." in constraint, it actually works! SCIP do not delete "obj_var" variable in transformed problem, and with createSol() and setSolVal() you can also set value of "obj_var" variable. But another thing, when you provide optimal solution SCIP runtime and number of proceeded nodes increase))) However at the output I can see that primal bound changed after providing optimal solution.

def create_TSP_model(c, V):
    model = Model("tsp")

    x = {}
    for i in V:
        for p in V:
            x[i,p] = model.addVar(vtype="B", name="x(%s,%s)"%(i,p))

    obj_var = model.addVar(vtype="C", name="obj_var")

    for i in V:
        model.addCons(quicksum(x[i, p] for p in V) == 1)
    for p in V:
        model.addCons(quicksum(x[i, p] for i in V) == 1)

    tmp = [i for i in V] + [0]
    model.addCons(obj_var + 1 >= quicksum(c[i, j]*quicksum(x[i, tmp[p]]*x[j, tmp[p+1]] for p in V) for i in V for j in V))

    model.setObjective(obj_var, "minimize")
    
    return model

EDIT: On problem with more variables it helps, just tried on TSP with 9 cities, and runtime without providing optimal solution 49 sec, and with 40 sec

EDIT2: But this "trick" with model actually give a bad performance compare to initial model where "obj_var" deleting. So, it is not the way to solve my problem:(

[Joao-Dionisio]

Hello, @SomovMike ! I am really sorry for taking so long to answer, I kept delaying and eventually, I forgot...

So it turns out that this is really a problem with SCIP, right? And your workaround didn't really solve the issue if I understood correctly.

[SomovMike]

Hello, unfortunately yes, I still can't solve this issue(

[Joao-Dionisio]

Did you manage to do something else, @SomovMike? What is the code you currently have?

[SomovMike]

Actually no, except one thing: SCIP made variable "obj var" multi aggregated, I don't know what it exactly means but you can set param in SCIP that it will not do that. And after that you can set value to variable "obj var" as well as other variables, and SCIP accepts such solution. But it works like a "trick" I described above, solution process become slower:(

[Joao-Dionisio]

When you say the solution process is slower, which process do you mean, @SomovMike?

From what I understand, you are solving a TSP instance to optimality, and then providing the optimal solution to a fresh copy of the model and then optimizing again, right? When you say it becomes slower, are you talking about the whole process, or just the second reoptimization?

If the entire process is slower, that is to be expected, as you are just repeating unnecessarily. However, just adding a solution to a model, is not a guarantee that it will be faster. One can construct scenarios where "artificially" introducing a solution makes it so other heuristics are called at a different time, or other methods don't work in the same way, and then the gain ends up being minimal or even nonexistent. Maybe this is what is happening.

[SomovMike]

No no no, I mean just Bnb process without proving any records. Turn this parameter off -- SCIP begins to solve the problem worse in general. But it is not a surprise, there is a reason why this parameter was originally enabled. But my goal is to provide records to SCIP without loosing in performance)

[Joao-Dionisio]

I don't think I understand what you are saying... What do you mean provide records in Bnb? Do you mean providing new solutions, updating the current best incumbent?

But you managed to solve the original problem, right, @SomovMike? Now it's more an issue related to your heuristic and the way SCIP deals with solutions provided by the user at runtime, yeah?

[SomovMike]

Sorry for misleading you. No, I didn't solve the initial issue. I still can't provide any solution to SCIP in a correct way.

[Joao-Dionisio]

No need to apologize, @SomovMike, I was the one who misunderstood. I also still can't solve your problem, sorry...

[Joao-Dionisio]

Hello, @SomovMike! What is the status of this issue? Shall we go for round 3 and see if we can figure out what's going on?

[Opt-Mucca]

So the requirement to fix this issue is probably: Add a Python function that wraps SCIPcreateOrigSol(), and add that to the exported SCIP functions.

[Joao-Dionisio]

@Opt-Mucca yeah, that seems like it would work, thanks!

But I think it would also be useful to have the possibility to transform an expression/solution/whatever from the original space into the transformed space and the other way around. What do you think?
This would need to be done in SCIP, something like SCIPgetTransformedExpression(original_expr) and SCIPgetOriginalExpressional(transformed_expr).