Below code is designed to reduce total supply chain costs, but it make production each month which has huge setup cost, instead it should hold inventory but it is not

can you help me in that?

import pandas as pd
import pulp
from pulp import LpVariable, LpProblem, lpSum, LpMinimize, LpInteger, LpBinary

filename = "Multi Period Multi Commodity SCNWD.xlsx"

# Read data from the file
demand_df = pd.read_excel(filename, sheet_name="Demand")
plant_capacity_df = pd.read_excel(filename, sheet_name="Plant Capacity")
warehouse_capacity_df = pd.read_excel(filename, sheet_name="Warehouse Capacity per Product")
fixed_plant_cost_df = pd.read_excel(filename, sheet_name="Fixed Plant Cost")
variable_production_cost_df = pd.read_excel(filename, sheet_name="Variable Production Cost")
inbound_cost_df = pd.read_excel(filename, sheet_name="Inbound Cost")
outbound_cost_df = pd.read_excel(filename, sheet_name="Outbound Cost")
fixed_warehouse_cost_df = pd.read_excel(filename, sheet_name="Fixed Warehouse Cost")
holding_cost_df = pd.read_excel(filename, sheet_name="Holding Cost")
setup_cost_df = pd.read_excel(filename, sheet_name="Setup Cost")

# Create dictionaries
demand = demand_df.set_index(["customer", "period", "product"])["demand"].to_dict()
plant_capacity = plant_capacity_df.set_index(["plant", "product"])["capacity"].to_dict()
warehouse_capacity = warehouse_capacity_df.set_index(["warehouse", "product"])["capacity"].to_dict()
fixed_plant_cost = fixed_plant_cost_df.set_index("plant")["cost"].to_dict()
variable_production_cost = variable_production_cost_df.set_index(["plant", "product"])["cost"].to_dict()
inbound_cost = inbound_cost_df.set_index(["plant","warehouse", "product"])["cost"].to_dict()
outbound_cost = outbound_cost_df.set_index(["warehouse",  "customer", "product"])["cost"].to_dict()
fixed_warehouse_cost = fixed_warehouse_cost_df.set_index("warehouse")["cost"].to_dict()
holding_cost = holding_cost_df.set_index("product")["cost"].to_dict()
setup_cost= setup_cost_df.set_index(["plant", "product"])["cost"].to_dict()

# Define Variables
customers= demand_df['customer'].unique().tolist()
plants= plant_capacity_df['plant'].unique().tolist()
warehouses= warehouse_capacity_df['warehouse'].unique().tolist()
products= demand_df['product'].unique().tolist()
periods= demand_df['period'].unique().tolist()

# Problem definition
problem = pulp.LpProblem("Multi-Commodity_Flow_Problem", pulp.LpMinimize)


# Decision variables
x = pulp.LpVariable.dicts("x", (plants, warehouses, products, periods), lowBound=0, cat="Continuous")
y = pulp.LpVariable.dicts("y", (warehouses, customers, products, periods), lowBound=0, cat="Continuous")
z_plant = pulp.LpVariable.dicts("z_plant", plants, cat="Binary")
z_warehouse = pulp.LpVariable.dicts("z_warehouse", warehouses, cat="Binary")
inventory = pulp.LpVariable.dicts("inv", (products, warehouses, range(0,13)), lowBound=0, cat="Continuous")

binary = pulp.LpVariable.dicts("binary", (plants, products, periods), lowBound=0, cat="Binary")
production_lot_size = pulp.LpVariable.dicts("lot_size", (plants, products, periods), lowBound=0, cat="Continuous")

# Objective
problem += (
    pulp.lpSum([fixed_plant_cost[p] * z_plant[p] for p in plants])
    + pulp.lpSum(
        [
            variable_production_cost[(p, pr)] * x[p][w][pr][t]
            for p in plants
            for w in warehouses
            for pr in products
            for t in periods
        ]
    )
    + pulp.lpSum(
        [
            inbound_cost[(p, w, pr)] * x[p][w][pr][t]
            for p in plants
            for w in warehouses
            for pr in products
            for t in periods
        ]
    )
    + pulp.lpSum(
        [
            outbound_cost[(w, c, pr)] * y[w][c][pr][t]
            for w in warehouses
            for c in customers
            for pr in products
            for t in periods
        ]
    )
    + pulp.lpSum([fixed_warehouse_cost[w] * z_warehouse[w] for w in warehouses])
    + pulp.lpSum(
        [
            binary[p][pr][t] * setup_cost[(p, pr)]
            for p in plants
            for pr in products
            for t in periods
        ]
    )
    + pulp.lpSum(
        [
            holding_cost[pr] * inventory[pr][w][t]
            for pr in products
            for t in periods
            for w in warehouses
        ]
   )
)



#----------------------- Constraints ----------------------
for p in plants:
    for pr in products:
        for t in periods:
            problem += production_lot_size[p][pr][t] >= binary[p][pr][t]

for pr in products:
    for t in periods:
        for w in warehouses:
            problem += (
                pulp.lpSum([x[p][w][pr][t] for p in plants])
                == pulp.lpSum([y[w][c][pr][t] for c in customers]) + inventory[pr][w][t]
            )

for c in customers:
    for pr in products:
        for t in periods:
            problem += (
                pulp.lpSum([y[w][c][pr][t] for w in warehouses])
                >= demand.get((c, t, pr), 0)
            )

for p in plants:
    for w in warehouses:
        for pr in products:
            for t in periods:
                problem += (
                    x[p][w][pr][t] <= plant_capacity.get((p, pr), 0) * z_plant[p]
                )

for w in warehouses:
    for pr in products:
        for t in periods:
            problem += (
                pulp.lpSum([y[w][c][pr][t] for c in customers])
                <= warehouse_capacity.get((w, pr), 0) * z_warehouse[w]
            )

# Additional constraints to link inventory and production_lot_size variables
for p in plants:
    for pr in products:
        for t in periods:
            problem += (
                pulp.lpSum([x[p][w][pr][t] for w in warehouses])
                == production_lot_size[p][pr][t]
            )

for pr in products:
    for w in warehouses:
        for t in range(1, len(periods)):
            problem += (
                inventory[pr][w][t - 1] + pulp.lpSum([production_lot_size[p][pr][t] for p in plants])
                - pulp.lpSum([demand.get((c, t, pr), 0) for c in customers])
                == inventory[pr][w][t]
            )

# Modified constraints to properly link inventory variable
for pr in products:
    for w in warehouses:
        for t in periods:
                # Start with initial inventory in the first period
                if t == min(periods):
                    problem += (
                        pulp.lpSum([x[p][w][pr][t] for p in plants for w in warehouses])
                        >= pulp.lpSum([y[w][c][pr][t] for w in warehouses for c in customers])
                        + inventory[pr][w][t]
                    )
                # For the rest of the periods, account for the previous period's inventory
                else:
                    problem += (
                        pulp.lpSum([x[p][w][pr][t] for p in plants for w in warehouses])
                        + inventory[pr][w][t - 1] - inventory[pr][w][t]
                        >= pulp.lpSum([demand.get((c, t, pr), 0) for c in customers])
                        
                    )


                    
# Solve the model with a time limit of 60 seconds
time_limit_in_seconds = 60 * 2

problem.solve(pulp.PULP_CBC_CMD(msg=1, timeLimit=time_limit_in_seconds))

# Print results
print("Objective value:", pulp.value(problem.objective))
print("Solution status:", pulp.LpStatus[problem.status])


# Print decision variables
for p in plants:
    for w in warehouses:
        for pr in products:
            for t in periods:
                if x[p][w][pr][t].varValue > 0:
                    print(f"x[{p}][{w}][{pr}][{t}] = {x[p][w][pr][t].varValue}")

for w in warehouses:
    for c in customers:
        for pr in products:
            for t in periods:
                if y[w][c][pr][t].varValue > 0:
                    print(f"y[{w}][{c}][{pr}][{t}] = {y[w][c][pr][t].varValue}")



for p in plants:
    if z_plant[p].varValue > 0:
        print(f"z_plant[{p}] = {z_plant[p].varValue}")

for w in warehouses:
    if z_warehouse[w].varValue > 0:
        print(f"z_warehouse[{w}] = {z_warehouse[w].varValue}")

# Print production lot size solution
print("Production lot size:")
for p in plants:
    for pr in products:
        for t in periods:
            if production_lot_size[p][pr][t].varValue > 0:
                print(f"Plant {p}, Product {pr}, Period {t}: {production_lot_size[p][pr][t].varValue}")

# Print inventory solution
for pr in products:
    for w in warehouses:
        for t in periods:
            if inventory[pr][w][t].varValue > 0:
                print(f"inv[{pr}][{t}] = {inventory[pr][w][t].varValue}")


        
#-------- Print the entire model output
for variable in problem.variables():
    if variable.varValue > 0:
        print(variable.name, "=", variable.varValue)


import matplotlib.pyplot as plt

def plot_product_results(product, time, demand, production, inventory):
    plt.figure(figsize=(12, 6))
    plt.title(f"Product: {product}")
    plt.plot(time, demand, marker='o', label="Demand")
    plt.bar(time, production, alpha=0.5, label="Production")
    plt.scatter(time, inventory, marker='s', color='red', label="Inventory")
    plt.xlabel("Month")
    plt.ylabel("Quantity")
    plt.legend()
    plt.grid()
    plt.xticks(time)
    plt.show()

    
    
# Initialize a dictionary to store the sum of demand for each product
sum_demand = {}

for pr in products:
    print(f"\n\nProduct: {pr}")
    print("Month".ljust(10) + "Production".ljust(15) + "Inventory".ljust(15) + "Total Demand")
    
    # Initialize lists to store demand, production, and inventory values for plotting
    demand_values = []
    production_values = []
    inventory_values = []
    
    # Initialize the sum of demand for the current product to 0
    sum_demand[pr] = 0
    
    for t in periods:
        # Calculate the total demand for all customers for the current product and time period
        total_demand = sum(demand[c, t, pr] for c in customers)
        
        # Add the total demand to the sum of demand for the current product
        sum_demand[pr] += total_demand
        
        # Initialize variables to store the production and inventory for the current product and time period
        total_production = 0
        total_inventory = 0
        
        # Calculate the total production and inventory for the current product and time period
        for p in plants:
            if production_lot_size[p][pr][t].varValue > 0:
                total_production += production_lot_size[p][pr][t].varValue
        for w in warehouses:
            if inventory[pr][w][t].varValue > 0:
                total_inventory += inventory[pr][w][t].varValue
        
        # Store the demand, production, and inventory values for plotting
        demand_values.append(total_demand)
        production_values.append(total_production)
        inventory_values.append(total_inventory)
        
        print(f"{t}".ljust(10) + f"{round(total_production,2)}".ljust(15) + f"{round(total_inventory, 2)}".ljust(15) + f"{total_demand}")
    
    # Plot the, demand, and inventory for the current product
    plot_product_results(pr, periods, demand_values, production_values, inventory_values)




 
# Define output dictionary
output = []

# Iterate over decision variables and add non-zero values to output dictionary
for p in plants:
    for w in warehouses:
        for pr in products:
            for t in periods:
                for c in customers:
                    
                        
                    if x[p][w][pr][t].varValue > 0:
                        # Add decision variable to output dictionary
                        output.append(
                            {
                                "plant": p,
                                "warehouse": w,
                                # "customer": c,
                                "product": pr,
                                "period": t,
                                "production": production_lot_size[p][pr][t].varValue,
                                "inbound_shipment": x[p][w][pr][t].varValue
                            }
                        )
                        
                    if y[w][c][pr][t].varValue > 0:
                        # Add decision variable to output dictionary
                        output.append(
                            {
                                "plant": p,
                                "warehouse": w,
                                "customer": c,
                                "product": pr,
                                "period": t,
                                "outbound_shipment": y[w][c][pr][t].varValue,
                            }
                        )

# Convert output to dataframe and export to Excel
output_df = pd.DataFrame(output)
print(output_df)