# -*- coding: utf-8 -*-
import pandas as pd
import pulp
from math import radians, sin, cos, sqrt, atan2

#  1. 读取数据
inventory_df = pd.read_excel('可调拨库存表.xlsx')
factory_loc_df = pd.read_excel('工厂位置.xlsx')
demand_df = pd.read_excel('需求表.xlsx')        # 型号、市、需求
city_loc_df = pd.read_excel('地名经纬度.xlsx')    # 市、经度、纬度

#  2. 数据合并
inventory_df['可调拨库存'] = pd.to_numeric(inventory_df['可调拨库存'], errors='coerce').fillna(0)
demand_df['需求'] = pd.to_numeric(demand_df['需求'], errors='coerce').fillna(0)

# 工厂库存合并
factory_capacity = {}
for _, row in inventory_df.groupby(['工厂', '型号'])['可调拨库存'].sum().reset_index().iterrows():
    f, p, cap = row['工厂'], row['型号'], int(row['可调拨库存'])
    if cap > 0:
        factory_capacity.setdefault(f, {})[p] = cap

# 城市需求合并
city_demand = {}
for _, row in demand_df.iterrows():
    city, p, dem = row['市'], row['型号'], int(row['需求'])
    if dem > 0:
        city_demand.setdefault(city, {})
        city_demand[city][p] = city_demand[city].get(p, 0) + dem

factory_lonlat = dict(zip(factory_loc_df['工厂'], zip(factory_loc_df['经度'], factory_loc_df['纬度'])))
city_lonlat    = dict(zip(city_loc_df['市'], zip(city_loc_df['经度'], city_loc_df['纬度'])))

# missing = set(city_demand) - set(city_lonlat)
# if missing:
#     raise ValueError(f"缺少经纬度：{missing}")

#  3. 距离
dist = {}
for f in factory_capacity:
    for c in city_demand:
        lon1, lat1 = factory_lonlat[f]
        lon2, lat2 = city_lonlat[c]
        dlon = radians(lon2-lon1)
        dlat = radians(lat2-lat1)
        a = sin(dlat/2)**2 + cos(radians(lat1))*cos(radians(lat2))*sin(dlon/2)**2
        dist[(f,c)] = 6371 * 2 * atan2(sqrt(a), sqrt(1-a))

#  4. 单位运费函数
def get_real_unit_cost(distance_km, product):
    d = distance_km
    p = str(product).strip()
    if p in ['GWBD-A2', 'GWBD-A3', 'GWBD-B']:
        if distance_km <= 100:
            return 15802.31
        elif distance_km <= 200:
            return 18445.12
        elif distance_km <= 300:
            return 22391.35
        elif distance_km <= 600:
            return 62.29 * d
        elif distance_km <= 900:
            return 44.75 * d
        elif distance_km <= 1200:
            return 38.23 * d
        elif distance_km <= 1500:
            return 36.08 * d
        elif distance_km <= 2000:
            return 33.44 * d
        elif distance_km <= 2500:
            return 32.76 * d
        elif distance_km <= 3500:
            return 26.68 * d
        else:
            return 24.27
    # 其他型号用第二套
    else:
        if distance_km <= 100:
            return 12859.46
        elif distance_km <= 200:
            return 14946.73
        elif distance_km <= 300:
            return 18662.46
        elif distance_km <= 600:
            return 52.83 * d
        elif distance_km <= 900:
            return 38.36 * d
        elif distance_km <= 1200:
            return 32.89 * d
        elif distance_km <= 1500:
            return 30.69 * d
        elif distance_km <= 2000:
            return 28.49 * d
        elif distance_km <= 2500:
            return 27.76 * d
        elif distance_km <= 3500:
            return 22.87 * d
        else:
            return 21.02 * d

#  5. 规划模型
prob = pulp.LpProblem("最小化运费", pulp.LpMinimize)

x = {}
for f in factory_capacity:
    for p in factory_capacity[f]:
        for city in city_demand:
            if city_demand[city].get(p, 0) > 0:
                x[(f,city,p)] = pulp.LpVariable(f"ship_{f}_{city}_{p}", lowBound=0, cat="Integer")

# 缺货惩罚
shortage = {(city,p): pulp.LpVariable(f"short_{city}_{p}", lowBound=0, cat="Continuous")
            for city in city_demand for p in city_demand[city]}

BIG_M = 1e10
prob += pulp.lpSum(
    x[k] * get_real_unit_cost(dist[(k[0],k[1])], k[2])   # 真实单位运费
    for k in x
) + pulp.lpSum(BIG_M * s for s in shortage.values())

# 库存硬约束
for f in factory_capacity:
    for p in factory_capacity[f]:
        prob += pulp.lpSum(x.get((f,c,p), 0) for c in city_demand) <= factory_capacity[f][p]

# 需求软约束
for city in city_demand:
    for p, dem in city_demand[city].items():
        prob += pulp.lpSum(x.get((f,city,p), 0) for f in factory_capacity) + shortage[(city,p)] >= dem

# 6. 求解
status = prob.solve(pulp.PULP_CBC_CMD(msg=True, timeLimit=600))

print(f"\n求解状态: {pulp.LpStatus[status]}")
if status != 1:
    print("无解（请检查总库存）")
else:
    print("有解")
    # real_cost = pulp.value(prob.objective)
    # print(f"实际运输总费用: {real_cost:,.2f} 元")

# ====================== 7. 输出 ======================
result = []
for (f,city,p), var in x.items():
    q = int(pulp.value(var) or 0)
    if q >= 1:
        d = dist[(f,city)]
        # unit = get_real_unit_cost(d, p)
        result.append({
            "工厂": f,
            "型号": p,
            "需求方": city,
            "调拨数量": q,
            # "距离_km": round(d, 1),
            # "真实单位运费": round(unit, 3),
            # "总运费": round(q * unit, 2)
        })

df = pd.DataFrame(result)
if not df.empty:
    df = df.sort_values(["工厂", "型号"])  # 按单位运费排序
    df[["工厂","型号","需求方","调拨数量"]].to_excel("调拨方案.xlsx", index=False)
    # print("\n调拨方案已保存（真实单位运费最优，优先满足单价更低的城市）")
    # print(df[["工厂","型号","需求方","调拨数量"]].to_string(index=False))
else:
    print("\n无调拨记录")

# 校验报告
# print("\n各城市需求满足情况：")
# for city in city_demand:
#     for p, dem in city_demand[city].items():
#         sup = sum(r["调拨数量"] for r in result if r["需求方"]==city and r["型号"]==p)
#         print(f"  {city} | {p}: 需求 {dem} → 调拨 {sup} → {'完美' if sup==dem else '异常'}")