I have written the following Python code for bus route optimisation. However, with any input parameters, the result is “Infeasible”. I think my mathematical formulation or code is not completely correct.
I considered minimization of walking distance to bus stations as the objective function. I also moved the bus route cost as a constraint and used linearization to replace the product of Sum(x_ij * sum(y_jk)). Can someone help me to fix the issue?
I have commented on every line, to make it clearer.

from pulp import LpMinimize, LpProblem, LpVariable, lpSum, LpBinary, LpStatus, value
import itertools

def network_design(Nodes, distance_matrix, cost, Budget):
    model = LpProblem("Transit_Network_Design", LpMinimize)
    N= len(Nodes) #we assume the first node is Depot
    Arcs = [(i, j) for i in range(N) for j in range(N)]
    # Decision variables
    x = LpVariable.dicts('x', Arcs, cat=LpBinary)  # If arc (i,j) is in transit route
    y = LpVariable.dicts('y', Arcs, cat=LpBinary)  # If demand centroid i is assigned to transit station j
    z = LpVariable.dicts('z', [(i, j, k) for i in range(N) for (j,k) in Arcs], lowBound = 0, upBound = 1)

    # Objective function
    model += lpSum(y[i, j] * distance_matrix[i][j] for i in range(N) for j in range(N) if (i,j) in Arcs)

    # Constraints
    model += lpSum(x[i, j] * cost[i][j] for i in range(N) for j in range(N) if (i, j) in Arcs) <= Budget  # Minimizing the total tour duration or cost

    # Auxiliary constraints for linearization
    for i in range(N):
        for j in range(N):
                for k in range(N):
                    if (j,k) in Arcs:
                      model += z[i, j, k] <= y[i, j]
                      model += z[i, j, k] <= x[j, k]
                      model += z[i, j, k] >= y[i, j] + x[j, k] - 1

    # Network flow constraints
    model += lpSum(x[i, j] for j in range(N) if i == 0) == 1  # A route starts from the depot
    model += lpSum(z[i, j, k] for i in range(N) for (j,k) in Arcs) == 1 # Each station must be served
    for i in range(N):
        model += lpSum(x[i, j] for j in range(N)) <= 1  # At most one arc can be in transit
        model += lpSum(x[i, j] for j in range(N)) - lpSum(x[j, i] for j in range(N)) == 0 
        
    model.solve()

    # Check if the problem was solved to optimality
    if LpStatus[model.status] != 'Optimal':
        return {
            "Status": LpStatus[model.status],
            "Message": "The problem did not solve to optimality. Check constraints and inputs."
        }

    # Prepare the results
    results = {
        "Status": LpStatus[model.status],
        "Objective Value": value(model.objective),
        "Selected Arcs": [(i, j) for (i, j) in Arcs if x[i, j].varValue is not None and x[i, j].varValue > 0.5],
        "Transit Station Connections": [(i, j) for i in range(N) for j in range(N) if y[i, j].varValue is not None and y[i, j].varValue > 0],
        "Selected z Variables": [(i, j, k) for i in range(N) for (j, k) in Arcs if z[i, j, k].varValue is not None and z[i, j, k].varValue > 0]
    }

    return results

# Example:
Nodes = ['Depot','A','B','C','D','E','F']
distance_matrix = [
    [1 , 29, 28, 16, 7 , 10, 15],
    [29, 1 , 17, 31, 22, 19, 14],
    [28, 17, 1 , 24, 21, 22, 17],
    [16, 31, 24, 1 , 9 , 13, 18],
    [7 , 22, 21, 9 , 1 , 3 , 8 ],
    [10, 19, 22, 13, 3 , 1 , 5 ],
    [15, 14, 17, 18, 8 , 5 , 1 ]
]
cost = [
    [1 , 10, 10, 10, 10, 10, 10],
    [10, 1 , 10, 10, 10, 10, 10],
    [10, 10, 1 , 10, 10, 10, 10],
    [10, 10, 10, 1 , 10, 10, 10],
    [10, 10, 10, 10, 1 , 10, 10],
    [10, 10, 10, 10, 10, 1 , 10],
    [10, 10, 10, 10, 10, 10, 1 ]
]

Budget = 100

results = network_design(Nodes, distance_matrix, cost, Budget)
for key, value in results.items():
    print(f"{key}: {value}")