Minimizing Cost Mathematical Linear Program

So I have this problem, and I have made a code in Gurobi and Python:

import gurobipy as gp
from gurobipy import GRB

# Define data
I = ['A', 'B', 'C']
M = list(range(1, 13))
S = ['Molding', 'Assembling', 'Finishing', 'Packaging']

# Demand in thousands of products
D = {
    ('A', 1): 10, ('B', 1): 11, ('C', 1): 10,
    ('A', 2): 20, ('B', 2): 13, ('C', 2): 20,
    ('A', 3): 25, ('B', 3): 12, ('C', 3): 50,
    ('A', 4): 30, ('B', 4): 11, ('C', 4): 10,
    ('A', 5): 50, ('B', 5): 10, ('C', 5): 20,
    ('A', 6): 70, ('B', 6): 10, ('C', 6): 50,
    ('A', 7): 10, ('B', 7): 10, ('C', 7): 10,
    ('A', 8): 12, ('B', 8): 12, ('C', 8): 20,
    ('A', 9): 10, ('B', 9): 11, ('C', 9): 50,
    ('A', 10): 300.2, ('B', 10): 10, ('C', 10): 10,
    ('A', 11): 20, ('B', 11): 50.4, ('C', 11): 20,
    ('A', 12): 14.8, ('B', 12): 10.6, ('C', 12): 55.5,
}

# Maximum production capacity per month for each production step
Cs = {
    'Molding': 80, 'Assembling': 100, 'Finishing': 100, 'Packaging': 250,
}

# Holding costs per 1000 products
H_finished = {'A': 800, 'B': 1500, 'C': 500}
H_packaged = {'A': 1000, 'B': 1000, 'C': 1000}

# Create a Gurobi model
model = gp.Model("Production_Schedule")

# Define decision variables with non-negativity constraints
X = {}
F = {}
R = {}
S_finished = {}
S_packaged = {}

for p in I:
    for m in M:
        for s in S:
            X[p, m, s] = model.addVar(vtype=GRB.CONTINUOUS, lb=0.0, name=f'X_{p}_{m}_{s}')
        F[p, m] = model.addVar(vtype=GRB.CONTINUOUS, lb=0.0, name=f'F_{p}_{m}')
        R[p, m] = model.addVar(vtype=GRB.CONTINUOUS, lb=0.0, name=f'R_{p}_{m}')
        S_finished[p, m] = model.addVar(vtype=GRB.CONTINUOUS, lb=0.0, name=f'S_{p}_finished_{m}')
        S_packaged[p, m] = model.addVar(vtype=GRB.CONTINUOUS, lb=0.0, name=f'S_{p}_packaged_{m}')

# Set objective function (minimize total holding costs)
model.setObjective(gp.quicksum(H_finished[p] * F[p, m] + H_packaged[p] * R[p, m] for p in I for m in M), GRB.MINIMIZE)

# Constraints
# 1. Demand Satisfaction
for p in I:
    for m in M:
        model.addConstr(gp.quicksum(X[p, m, s] for s in S) + R[p, m] == D[p, m], f'Demand_Satisfaction_{p}_{m}')

# 2. Production Capacity
for s in S:
    model.addConstr(gp.quicksum(X[p, m, s] for p in I for m in M) <= Cs[s], f'Production_Capacity_{s}')

# 3. Inventory Balance
# Inventory Balance for m = 1
for p in I:
    model.addConstr(F[p, 1] == gp.quicksum(X[p, 1, s] for s in S) - D[p, 1], f'Inventory_Balance_Finished_{p}_1')
    model.addConstr(R[p, 1] == 0, f'Inventory_Balance_Packaged_{p}_1')  # Assuming no packaged inventory at the start

# Inventory Balance for m > 1
for p in I:
    for m in M[2:]:
        model.addConstr(F[p, m] == F[p, m - 1] + gp.quicksum(X[p, m - 1, s] for s in S) - D[p, m], f'Inventory_Balance_Finished_{p}_{m}')
        model.addConstr(R[p, m] == S_packaged[p, m - 1], f'Inventory_Balance_Packaged_{p}_{m}')

# 4. Production Order
for p in I:
    for m in M:
        model.addConstr(X[p, m, 'Molding'] <= X[p, m, 'Assembling'], f'Production_Order_Molding_Assembling_{p}_{m}')
        model.addConstr(X[p, m, 'Assembling'] <= X[p, m, 'Finishing'], f'Production_Order_Assembling_Finishing_{p}_{m}')
        model.addConstr(X[p, m, 'Finishing'] <= X[p, m, 'Packaging'], f'Production_Order_Finishing_Packaging_{p}_{m}')

# 5. First 3 Production Steps
for p in I:
    for m in M:
        molding_assembling_finishing_sum = gp.quicksum(X[p, m, s] for s in S[:3])
        individual_step_sum = X[p, m, 'Molding'] + X[p, m, 'Assembling'] + X[p, m, 'Finishing']
        model.addConstr(molding_assembling_finishing_sum == individual_step_sum, f'First_3_Production_Steps_{p}_{m}')

# 6. Storage of Finished Product
for p in I:
    for m in M[:-1]:
        model.addConstr(S_finished[p, m] == 0, f'Storage_Finished_Product_{p}_{m}')

# 7. Storage for Demand
for p in I:
    for m in M:
        model.addConstr(S_packaged[p, m] <= D[p, m], f'Storage_For_Demand_{p}_{m}')

# Set MIPGap parameter
model.setParam('MIPGap', 0)

# Solve the model
model.optimize()

# Check the status and retrieve the optimal solution if available
if model.status == GRB.OPTIMAL:
    # Optimal cost value
    optimal_cost = model.objVal
    print(f"Optimal Cost: {optimal_cost}")

    # Optimal production schedule (values of X)
    production_schedule = {(p, m, s): X[p, m, s].x for p in I for m in M for s in S}
    print("Optimal Production Schedule (X):")
    for (p, m, s), value in production_schedule.items():
        if value > 0:
            print(f"X_{p}_{m}_{s} = {value}")

    # You can similarly retrieve the values of other decision variables like F, R, S_finished, and S_packaged
else:
    print("No optimal solution found.")

# Dispose of the model
model.dispose()

But the output still says no optimal solution found? Can anyone give me some feedback on this? Thank you. I also would like the output to be in a matrix like month and the steps