Search code examples
pythonlinear-programminggekkointeger-programming

setting a lower bound constraint based on condition in gekko


I am trying to set up something similar to the facility location problem using gekko in python. I'm trying to minimize the distance between a set of facilities and corresponding counties while ensuring that if a facility is used as an optimal assignment it accounts for at least 10% of the total population in all of the counties.

I'm trying to allow the algorithm to decide the best facilities to use, so I dont want to limit how many it can select.

so, my objective is to minimize the distance between the chosen facilities and the counties that are mapped to it. my constraints are that each county can only be mapped to one facility. And, if a facility is selected in the final solution, the sum of the population in each county mapped to that facility must be at least 10% of the total population in all the counties.

I keep getting "@error: Solution Not Found" and dont know other steps to take using this method. So, is there a better way to set this problem up, and if not what other algorithms would you recommend using to solve this problem?

Relevant Data

  • Note: this is a very small selection of my real dataset, I'm just trying to show what it looks like

  • county_idx = list of integers from 0 to I-1 where I = total # of counties

  • selected_facilities = list of integers from 0 to J-1 where J = total # of facilities

  • distances_mapped = array where rows are the counties index and columns are the facility index. The values are corresponding distances between the county and facility

  • populations_mapped = pandas series that has the population of each county where the rows are indexed by the county_idx and value is the corresponding population


county_idx = [0,1,2]
selected_facilities = [0, 1, 2]
distances_mapped = np.array([[193, 85, 226],
                    [139, 112, 241],
                    [175, 110, 249]])

populations_mapped = [981447, 327286, 176622]

cutoff_population = 153

Relevant Code

m = GEKKO(remote=False)
n_rows = len(county_idx)
n_cols = len(selected_facilities)

# Create binary variables for facility selection

facility_assignment = m.Array(m.Var, (n_cols), lb=0, ub=1, integer=True)

# Create binary variables for county-facility mapping
assignment_matrix = m.Array(m.Var, (n_rows, n_cols), lb=0, ub=1, integer=True)

# Define the objective function - weighted sum of distances
m.Minimize(m.sum(assignment_matrix * distances_mapped))

# Each county is serviced by exactly one facility
for i in county_idx:
    m.Equation(m.sum([assignment_matrix[(i, j)] for j in selected_facilities]) == 1)

# need to bind facility assignment to make sure it ties to the assignment matrix
for j in selected_facilities:
    # need to make sure facility assignment updates when a facility is selected
    # so, if any of the counties i have a facility j mapped to it then facility_assignement needs to keep track of that
    m.Equation(facility_assignment[j] == m.sum([assignment_matrix[(i, j)] for i in county_idx]))


# Ensure that the population assigned to a facility is >= x% of the total population
for j in selected_facilities:
    # logic: if a facility is selected then the sum of its population must be >= x% of the total population
    m.Equation(m.sum([assignment_matrix[(i, j)] * populations_mapped.iloc[i, 0] for i in county_idx]) >= cutoff_population * facility_assignment[j])

# ensures integer solutions using apopt solver
m.options.SOLVER = 1

m.options.IMODE = 3

# Solve the optimization problem
m.solve(disp=True, debug=True)

Solution

  • You are very close. You need to change your second constraint to a "Big M" constraint. Right now, you have this:

        m.Equation(facility_assignment[j] == m.sum([assignment_matrix[(i, j)] for i in county_idx]))
    

    which has the unintended consequence of allowing the sum of the matrix for the selected facility to be 1 or 0, meaning it is limiting the assignment to a single county.

    What you want is to make the facility_assignment variable 1 if one or more assignments occur, so that is the job for "big M", in this case, M would logically be the number of counties

    M = len(county_idx)
    

    and in the constraint:

        m.Equation(facility_assignment[j] * M >= m.sum([assignment_matrix[(i, j)] for i in county_idx]))
    

    which will force facility_assigned to be 1 if 1, 2, 3, ... M county assignments are given to it, and trigger the enforcement of your population constraint, which is correctly written (assuming you define cutoff population correctly)

    As a nit-pick, I would also rename facility_assigned to facility_used because I think it says more clearly that it is used (at least once) and the assignment piece is handled in your other variable.