Here is the homework prompt:

Try to satisfy population balance (one-person-one-vote). That is, congressional districts should have approximately the same population. Consider strategies for assigning more than one representative to counties with high-population centers as long as elections are county-wide. Do not divide counties geographically.
Use integer programming (set partitioning) to obtain an algorithmic/optimal redistricting. Assign every county in your selected state to exactly one congressional district while striving to meet your objective through maximization or minimization.
Try to ensure that congressional districts are composed of counties that are adjacent to one another. Describe constraints or objectives employed to accomplish this goal. Note any difficulties encountered in setting up constraints or objectives.
Try to ensure that congressional districts are geographically compact.
Solve the integer programming problem using Python PuLP.  Note any difficulties encountered, given the size of the integer programming problem.
Consider secondary goals of redistricting, such as encouraging equal representation across races. For example, you may try to achieve as much racial balance (percentage white alone versus other races) as possible across all congressional districts.
 Draw color-coded maps for your algorithmic/optimal redistricting and for the actual redistricting that was implemented by your selected state
I have attached my table and code. I am getting an infeasible solution, how do I solve that? Do I need to add more constraints? Additionally, I do not know how I would map the redistricting solution in python. Please help!

countypopulationlonglatFIPS
0Adams County20961-118.56035646.98352453001
1Asotin County22508-117.20279146.19177153003
2Benton County212791-119.51133446.23976053005
3Chelan County79926-120.61882847.86921953007
4Clallam County77805-123.92771148.04911753009
5Clark County516779-122.48231145.77919853011
6Columbia County4026-117.90782546.29796153013
7Cowlitz County111956-122.68027346.19321553015
8Douglas County44192-119.69170747.73589653017
9Ferry County7448-118.51632448.47039753019
10Franklin County98678-118.89847546.53504653021
11Garfield County2363-117.54464846.43131053023
12Grant County101311-119.45139647.20582653025
13Grays Harbor County77038-123.77304847.14971053027
14Island County86625-122.54939248.16376153029
15Jefferson County33589-123.59359247.74901353031
16King County2266789-121.80636147.49034853033
17Kitsap County277673-122.67340247.61342453035
18Kittitas County45189-120.67971247.12448653037
19Klickitat County23271-120.78950345.87361553039
20Lewis County85370-122.39334446.57770753041
21Lincoln County11601-118.41865647.57605453043
22Mason County68166-123.19203447.34793253045
23Okanogan County43127-119.74060348.54863053047
24Pacific County24113-123.70572346.55548053049
25Pend Oreille County14179-117.27383248.53235653051
26Pierce County927380-122.11221247.02680653053
27San Juan County18662-122.95346348.57856353055
28Skagit County131179-121.73192348.47967853057
29Skamania County12460-121.91489146.02281053059
30Snohomish County840079-121.69776048.04772453061
31Spokane County549690-117.40394147.62078053063
32Stevens County48229-117.85507148.39899253065
33Thurston County298758-122.83306546.92672953067
34Wahkiakum County4688-123.42476046.29242253069
35Walla Walla County61890-118.47853346.22968353071
36Whatcom County230677-121.72095248.82595353073
37Whitman County47619-117.52297746.90119153075
38Yakima County257001-120.73873946.45687553077
 

import pandas as pd
import pulp
import numpy as np
from tqdm import tqdm

df = pd.read_excel("counties-table.xlsx")
df = df[['county', 'population', 'long', 'lat' ,'FIPS']]

df

import math
from math import pi, pow, sin, cos, asin, sqrt, floor
from scipy import stats
from pyproj import Proj


def degrees_to_radians(x):
     return((np.pi/180)*x)
     
def lon_lat_distance_miles(lon_a,lat_a,lon_b,lat_b):
    radius_of_earth = 24872/(2*np.pi)
    c = np.sin((degrees_to_radians(lat_a) - \
    degrees_to_radians(lat_b))/2)**2 + \
    np.cos(degrees_to_radians(lat_a)) * \
    np.cos(degrees_to_radians(lat_b)) * \
    np.sin((degrees_to_radians(lon_a) - \
    degrees_to_radians(lon_b))/2)**2
    return(2 * radius_of_earth * (np.arcsin(np.sqrt(c))))    

def lon_lat_distance_meters (lon_a,lat_a,lon_b,lat_b):
    return(lon_lat_distance_miles(lon_a,lat_a,lon_b,lat_b) * 1609.34)


counties_Lat_Dic = dict(zip(df.county, df.lat))
counties_Long_Dic = dict(zip(df.county, df.long))
counties_Pop_Dic = dict(zip(df.county, df.population))

#creating a data frame with lat and long values and pairwise distance between counties
latlon = df[['lat', 'long']].values
n = df.shape[0]
dists = np.zeros((n,n))
for i in range(n):
    for j in range(i+1,n):
        dists[i,j] = lon_lat_distance_meters(latlon[i,1], latlon[i,0], latlon[j,1], latlon[j,0])
        dists[j,i] = dists[i,j]
       
dists_df = pd.DataFrame(dists, columns=df['county'].values, index=df['county'])
pops_df = df.set_index('county')['population'].astype(float)

dists_df

max_distance = 200000
max_district = 10
max_district_size = 5
max_pop_diff = 402488  #standard deviation of population

counties = df['county'].values.tolist()[:39]

#Define compactness function
def compactness(district):

    k = len(district)
    if k == 1:
        return 0
    cur = dists_df.loc[district, district].values
    return cur[np.triu_indices(k, 1)].mean()

def distance(district):
    return np.max(dists_df.loc[district, district].values) <= max_distance

def population(district):
    min_pop = min(pops_df.loc[list(district)])
    max_pop = max(pops_df.loc[list(district)])
    return max_pop - min_pop <= max_pop_diff

#create list of all possible districts
possible_districts = []
for county in tqdm(counties):
    neighbors = dists_df.loc[county,counties]
    neighbors = neighbors[neighbors <= max_distance].index.tolist()
    possible_districts += [tuple(sorted(c)) for c in pulp.allcombinations(neighbors, max_district_size) if distance(c) and population(c)]
possible_districts = list(set(possible_districts))
print(len(possible_districts))

#create a binary variable to state that a district is used
x = pulp.LpVariable.dicts('district', possible_districts,
                            lowBound = 0,
                            upBound = 1,
                            cat = pulp.LpInteger)

redistrict_model = pulp.LpProblem("Redistricting Model", pulp.LpMinimize)

redistrict_model += sum([compactness(district) * x[district] for district in possible_districts])

#specify the maximum number of districts
redistrict_model += sum([x[district] for district in possible_districts]) <= max_district, "Maximum_number_of_districts"


#A county can be assigned to one and only one district
for county in counties:
    redistrict_model += sum([x[district] for district in possible_districts
                                if county in district]) == 1, "Must_zone_%s"%county

   
redistrict_model.solve()

if redistrict_model.status != -1:
    print("The choosen districts are out of a total of %s:"%len(possible_districts))
    for district in possible_districts:
        if x[district].value() == 1.0:
            print(district)