Pyomo: Cant transform non linear constraint to a linear one

Question

I am currently working on a linear programming problem,where we assign sessions to certain caregivers, and a constraint gives me an error for non linearity. Before I put the code for the constraint, I will introduce all the model items used in such constraint.

• model.CASE_COMBINATIONS is a set which corresponds to the cross product between two series with the ids of the sessions(also a set). This is supposed the represent all combinations of cases. Here is the code for this set:

    model.CASE_COMBINATIONS = pe.Set(
            initialize=product(
                self.df_sessions["idx"].unique(),
                self.df_sessions["idx"].unique(),
            ),
            dimen=2,
        )

• model.SESSION_ASSIGNED is a variable of the model, which indicates if a certain session is selected to be performed by a certain caregiver. This variables does this by assigning 1 to the tuple(session, caregiver) if the caregiver is assigned to that session, and 0 otherwise.

Here is the code for model.SESSION_ASSIGNED and model.TASKS:

model.SESSION_ASSIGNED = pe.Var(model.TASKS, domain=pe.Binary)

model.TASKS = pe.Set(
            initialize=model.CASES * model.CAREGIVERS, dimen=2
        )

• model.CASE_BIGGER is a parameter which gives us, for each tuple of sessions(case1, case2) if the idx of case1 is bigger than case2. We do this because the idx are ordered, so this allows us to know which session occurs first. Here is the code:

    model.CASE_BIGGER = pe.Param(
                model.CASE_COMBINATIONS, initialize=self._generate_case_bigger()

    def _generate_case_bigger(self):
            case_bigger = {}
            for case1, case2 in product(
                self.df_sessions["idx"].unique(), self.df_sessions["idx"].unique()
            ):
                case_bigger[(case1, case2)] = int(case1 > case2)
            return case_bigger

• model.COMMUTE corresponds to a parameter with the commute time between two different locations of clients, where for each tuple (session1, session2), we will have a mapping to the commute time between the two. Here is the code:

    model.COMMUTE = pe.Param(
                model.CLIENT_CONNECTIONS,
                initialize=self._generate_clients_commute(),
            )

• model.IDX_CLIENTS just corresponds is just a mapping between two different sets of ids that we have for sessions, and I believe it is mostly irrelevant to the question at hand.

Now for the constraint with non linearity issues:

def commute_care(model, caregiver):
    commute_expr = sum(
        [
            model.SESSION_ASSIGNED[case[0], caregiver] *
            model.SESSION_ASSIGNED[case[1], caregiver] *
            ( 
                model.CASE_BIGGER[case[1], case[0]] *
                model.COMMUTE[(model.IDX_CLIENTS[case[0]], model.IDX_CLIENTS[case[1]])] +

                (1 - model.CASE_BIGGER[case[1], case[0]]) *
                model.COMMUTE[(model.IDX_CLIENTS[case[1]], model.IDX_CLIENTS[case[0]])]
            )
            for case in model.CASE_COMBINATIONS
        ]
    )
    return model.COMMUTE_CARE[caregiver] == commute_expr

With this constraint, we want to calculate the commute time for each caregiver. We iterate over the cross product of the different cases, and we calculate, for each pair of cases, if they were assigned to the same caregiver and if so, we calculate the commute time between them. The sum of these commute times is then equal to the sum of all these commute times. Does anyone have any idea on how to transform this constraint into a linear constraint? Thank you

Answer 1

Taking a step back from the code, you have expressions of the form:

a[i,j]*x[i]*x[j]

where a is constant, and x are binary variables. This can be linearized using a well-known reformulation:

a[i,j]*y[i,j]
y[i,j] <= x[i]
y[i,j] <= x[j]
y[i,j] >= x[i]+x[j]-1

where y is an extra binary variable. There are some extra things we can exploit:

• y can be relaxed to be continuous between 0 and 1
• depending on the objective we may drop some of the new constraints as they will automatically hold
• we can exploit symmetry between i,j