import pandas as pdimport pulp as plFILE = "lp_realistic_data_v2.xlsx"sites_df   = pd.read_excel(FILE, sheet_name="sites")mkts_df    = pd.read_excel(FILE, sheet_name="markets")prod_df    = pd.read_excel(FILE, sheet_name="products")demand_df  = pd.read_excel(FILE, sheet_name="demand")ship_df    = pd.read_excel(FILE, sheet_name="shipping")bom_df     = pd.read_excel(FILE, sheet_name="bom")rlim_df    = pd.read_excel(FILE, sheet_name="resource_limits")upg_df     = pd.read_excel(FILE, sheet_name="site_upgrades")params_df  = pd.read_excel(FILE, sheet_name="params")sites    = sites_df["site"].tolist()markets  = mkts_df["market"].tolist()products = prod_df["product"].tolist()fixed_cost = dict(zip(sites_df["site"], sites_df["fixed_cost"]))base_cap   = dict(zip(sites_df["site"], sites_df["base_capacity"]))site_addon = dict(zip(sites_df["site"], sites_df["prod_cost_addon"]))mkt_weight = dict(zip(mkts_df["market"], mkts_df["weight"]))price          = dict(zip(prod_df["product"], prod_df["unit_price"]))base_unit_cost = dict(zip(prod_df["product"], prod_df["base_unit_cost"]))demand    = {(r["market"], r["product"]): float(r["demand"]) for _, r in demand_df.iterrows()}ship_cost = {(r["site"], r["market"]): float(r["ship_cost"]) for _, r in ship_df.iterrows()}labor_use = dict(zip(bom_df["product"], bom_df["labor_per_unit"]))mat_use   = dict(zip(bom_df["product"], bom_df["material_per_unit"]))res_limit = dict(zip(rlim_df["resource"], rlim_df["limit"]))params = dict(zip(params_df["param"], params_df["value"]))K              = int(params["K"])INVEST_BUDGET  = float(params["INVEST_BUDGET"])SERVICE_LEVEL  = float(params["SERVICE_LEVEL"])ALLOW_BACKORDER = int(params["ALLOW_BACKORDER"])SHORT_PENALTY  = float(params.get("SHORT_PENALTY", 0.0))TOTAL_DEMAND = sum(demand[(m, p)] for m in markets for p in products)mkt_total_demand = {m: sum(demand[(m, p)] for p in products) for m in markets}

m1 = pl.LpProblem("Location", pl.LpMinimize)open_site = pl.LpVariable.dicts("open", sites, 0, 1, cat=pl.LpBinary)m1 += (    pl.lpSum(fixed_cost[s] * open_site[s] for s in sites)    + pl.lpSum(        mkt_weight[m] * mkt_total_demand[m] * ship_cost[(s, m)] * open_site[s]        for s in sites for m in markets    ))m1 += pl.lpSum(open_site[s] for s in sites) == Km1.solve(pl.PULP_CBC_CMD(msg=False))selected_sites = [s for s in sites if pl.value(open_site[s]) > 0.5]print("Status:", pl.LpStatus[m1.status])print("Selected sites:", selected_sites)print("Objective (location):", pl.value(m1.objective))print()if pl.LpStatus[m1.status] != "Optimal":    raise RuntimeError("Part 1 did not solve to optimal. Check data/constraints.")

m2 = pl.LpProblem("Investment", pl.LpMaximize)upg_pairs = [(r["site"], r["upgrade"]) for _, r in upg_df.iterrows()]upg_cost  = {(r["site"], r["upgrade"]): float(r["upgrade_cost"]) for _, r in upg_df.iterrows()}upg_gain  = {(r["site"], r["upgrade"]): float(r["capacity_gain"]) for _, r in upg_df.iterrows()}x = pl.LpVariable.dicts("upgrade", upg_pairs, 0, 1, cat=pl.LpBinary)m2 += pl.lpSum(upg_gain[pair] * x[pair] for pair in upg_pairs)for (s, u) in upg_pairs:    m2 += x[(s, u)] <= (1 if s in selected_sites else 0)for s in selected_sites:    m2 += pl.lpSum(x[(ss, u)] for (ss, u) in upg_pairs if ss == s) <= 1m2 += pl.lpSum(upg_cost[pair] * x[pair] for pair in upg_pairs) <= INVEST_BUDGETm2.solve(pl.PULP_CBC_CMD(msg=False))chosen_upgrades = [pair for pair in upg_pairs if pl.value(x[pair]) > 0.5]final_cap = {}for s in sites:    cap = base_cap[s] * (1 if s in selected_sites else 0)    for (ss, u) in chosen_upgrades:        if ss == s:            cap += upg_gain[(ss, u)]    final_cap[s] = capTOTAL_CAP = sum(final_cap[s] for s in selected_sites)INV_USED = sum(upg_cost[pair] for pair in chosen_upgrades)print("Status:", pl.LpStatus[m2.status])print("Chosen upgrades:", chosen_upgrades if chosen_upgrades else "None")print(f"Invest used: {INV_USED} / {INVEST_BUDGET}")print("Final capacity (selected sites):", {s: final_cap[s] for s in selected_sites})print("Total final capacity:", TOTAL_CAP)print()if pl.LpStatus[m2.status] != "Optimal":    raise RuntimeError("Part 2 did not solve to optimal. Check data/constraints.")

m3 = pl.LpProblem("Production_and_Shipping", pl.LpMaximize)ship = pl.LpVariable.dicts(    "ship",    [(s, m, p) for s in selected_sites for m in markets for p in products],    lowBound=0,    cat=pl.LpContinuous)short = pl.LpVariable.dicts(    "short",    [(m, p) for m in markets for p in products],    lowBound=0,    cat=pl.LpContinuous)m3 += pl.lpSum(    price[p] * ship[(s, m, p)]    - (base_unit_cost[p] + site_addon[s]) * ship[(s, m, p)]    - ship_cost[(s, m)] * ship[(s, m, p)]    for s in selected_sites for m in markets for p in products) - SHORT_PENALTY * pl.lpSum(short[(m, p)] for m in markets for p in products)for m in markets:    for p in products:        m3 += pl.lpSum(ship[(s, m, p)] for s in selected_sites) + short[(m, p)] == demand[(m, p)]if ALLOW_BACKORDER == 0:    for m in markets:        for p in products:            m3 += short[(m, p)] == 0for s in selected_sites:    m3 += pl.lpSum(ship[(s, m, p)] for m in markets for p in products) <= final_cap[s]m3 += pl.lpSum(labor_use[p] * ship[(s, m, p)] for s in selected_sites for m in markets for p in products) <= res_limit["labor"]m3 += pl.lpSum(mat_use[p]   * ship[(s, m, p)] for s in selected_sites for m in markets for p in products) <= res_limit["material"]m3 += pl.lpSum(ship[(s, m, p)] for s in selected_sites for m in markets for p in products) >= SERVICE_LEVEL * TOTAL_DEMANDm3.solve(pl.PULP_CBC_CMD(msg=False))