#!/usr/bin/env python3 """ inventory_transport_optimizer.py ================================= 成品油油库库存-运输协同优化求解器 实现内容: 1. 确定性等价模型 (CE) — 期望值近似 2. 场景法随机规划模型 (SP) — Scenario-Based DEP 3. 混合 GA-SA 启发式求解 — 大规模场景 4. MPC 滚动时域优化 — 在线运营 5. 基线方法对比 — (s,S), 固定批量, 经验规则 依赖: numpy, openpyxl (数据读写) 无商业求解器依赖，纯 Python 实现。 用法: python inventory_transport_optimizer.py """ import sys, os sys.path.insert(0, '/usr/lib/python3/dist-packages') import numpy as np from copy import deepcopy import time, warnings warnings.filterwarnings('ignore') # ═══════════════════════════════════════════════════════════ # 0. 数据加载 # ═══════════════════════════════════════════════════════════ def load_data_from_excel(filepath): """从 Excel 加载模型参数与需求数据""" from openpyxl import load_workbook wb = load_workbook(filepath, data_only=True) # 从参数配置Sheet读取 s1 = wb['参数配置'] params = {} for r in range(2, s1.max_row + 1): key = s1.cell(row=r, column=2).value val = s1.cell(row=r, column=3).value if key: params[key] = val # 补给源参数 s2 = wb['补给源参数'] sources = {} for r in range(2, s2.max_row + 1): sid = int(s2.cell(row=r, column=1).value) sources[sid] = { 'name': s2.cell(row=r, column=2).value, 'Q_max': float(s2.cell(row=r, column=3).value), 'L_km': float(s2.cell(row=r, column=4).value), 'Ct': float(s2.cell(row=r, column=5).value), } # 需求数据 T=30 s3 = wb['需求_T30'] demand_30 = np.array([float(s3.cell(row=r, column=2).value) for r in range(2, s3.max_row+1)]) # 场景数据 s6 = wb['场景_T30_前10'] n_scenes = s6.max_column - 1 T = s6.max_row - 1 scenarios = np.zeros((T, n_scenes)) for t in range(T): for s in range(n_scenes): scenarios[t, s] = float(s6.cell(row=t+2, column=s+2).value) # 聚合场景 try: s7 = wb['聚合场景_K30'] K = s7.max_column - 1 T2 = s7.max_row - 2 clust = np.zeros((T2, K)) probs = np.zeros(K) for k in range(K): probs[k] = float(s7.cell(row=2, column=k+2).value) for t in range(T2): for k in range(K): clust[t, k] = float(s7.cell(row=t+3, column=k+2).value) except: clust, probs = None, None return { 'M': int(params.get('M', 2)), 'T': int(params.get('T_s', 30)), 'V': float(params.get('V', 12000)), 'I0': float(params.get('I0', 3500)), 'theta': float(params.get('θ', 0.003)), 'alpha': float(params.get('α', 0.95)), 'Ch': float(params.get('Ch', 12)), 'Cs': float(params.get('Cs', 850)), 'Cl': float(params.get('Cl', 7200)), 'sources': sources, 'demand': demand_30, 'scenarios': scenarios, 'clustered': clust, 'cluster_probs': probs, } # ═══════════════════════════════════════════════════════════ # 1. 模型核心：成本计算 # ═══════════════════════════════════════════════════════════ class InventoryTransportModel: """库存-运输协同优化模型""" def __init__(self, data): self.M = data['M'] self.T = data['T'] self.V = data['V'] self.I0 = data['I0'] self.theta = data['theta'] self.alpha = data['alpha'] self.Ch = data['Ch'] self.Cs = data['Cs'] self.Cl = data['Cl'] self.sources = data['sources'] self.Q_max = np.array([data['sources'][i]['Q_max'] for i in range(self.M)]) self.L = np.array([data['sources'][i]['L_km'] for i in range(self.M)]) self.Ct = np.array([data['sources'][i]['Ct'] for i in range(self.M)]) def compute_cost_deterministic(self, Q, demand): """ 确定性场景下计算总成本 Q: (M, T) 运输批量矩阵 demand: (T,) 需求序列 返回: (total_cost, holding, shortage, transport, loss) """ M, T = Q.shape I = np.zeros(T) holding = 0; shortage = 0; transport = 0; loss = 0 I_prev = self.I0 for t in range(T): arrival = np.sum(Q[:, t]) * (1 - self.theta) S_t = I_prev + arrival served = min(S_t, demand[t]) I[t] = max(0, S_t - demand[t]) # 各项成本 holding += self.Ch * I[t] shortage += self.Cs * max(0, demand[t] - S_t) for i in range(M): transport += self.Ct[i] * self.L[i] * Q[i, t] loss += self.Cl * self.theta * Q[i, t] # 库容约束检查 if I[t] > self.V: holding += 1e6 * (I[t] - self.V) # 强惩罚 I_prev = I[t] total = holding + shortage + transport + loss return total, holding, shortage, transport, loss def compute_cost_scenario(self, Q, scenarios, probs=None): """ 场景法下计算期望总成本 Q: (M, T) 运输批量矩阵 scenarios: (T, S) 需求场景矩阵 probs: (S,) 场景概率 """ M, T = Q.shape _, S = scenarios.shape if probs is None: probs = np.ones(S) / S total_exp = 0; hold_exp = 0; short_exp = 0 trans = 0; loss = 0 # 运输成本与损耗与场景无关 for t in range(T): for i in range(M): trans += self.Ct[i] * self.L[i] * Q[i, t] loss += self.Cl * self.theta * Q[i, t] # 按场景计算库存与缺货 for s in range(S): I = np.zeros(T) I_prev = self.I0 for t in range(T): arrival = np.sum(Q[:, t]) * (1 - self.theta) S_t = I_prev + arrival I[t] = max(0, S_t - scenarios[t, s]) hold_exp += probs[s] * self.Ch * I[t] short_exp += probs[s] * self.Cs * max(0, scenarios[t, s] - S_t) I_prev = I[t] total = hold_exp + short_exp + trans + loss return total, hold_exp, short_exp, trans, loss def service_level(self, Q, scenarios, probs=None): """计算需求满足率""" M, T = Q.shape _, S = scenarios.shape if probs is None: probs = np.ones(S) / S total_demand = 0; total_shortage = 0 for s in range(S): I_prev = self.I0 for t in range(T): arrival = np.sum(Q[:, t]) * (1 - self.theta) S_t = I_prev + arrival served = min(S_t, scenarios[t, s]) total_demand += probs[s] * scenarios[t, s] total_shortage += probs[s] * max(0, scenarios[t, s] - S_t) I_prev = max(0, S_t - scenarios[t, s]) return 1.0 - total_shortage / max(total_demand, 1e-6) # ═══════════════════════════════════════════════════════════ # 2. 基线方法 # ═══════════════════════════════════════════════════════════ def baseline_ss_policy(model, demand_mu, demand_std): """(s, S) 策略""" M, T = model.M, model.T Q = np.zeros((M, T)) I_prev = model.I0 # 计算 (s, S) 参数 z = 1.645 # α=0.95 avg_daily = demand_mu ss = z * demand_std * np.sqrt(2) # 安全库存 s = avg_daily * 3 + ss # 订货点 S_level = avg_daily * 5 + ss # 订货上限 for t in range(T): arrival = np.sum(Q[:, t]) * (1 - model.theta) I_after = I_prev + arrival - (demand_mu if t < len(np.zeros(T)) else 0) if I_after < s: order_qty = min(S_level - I_after, model.Q_max[0]) Q[0, t] = order_qty I_prev = max(0, I_after) return Q def baseline_fixed_batch(model, batch_size=200, interval=3): """固定批量补货策略""" M, T = model.M, model.T Q = np.zeros((M, T)) for t in range(0, T, interval): Q[0, t] = min(batch_size, model.Q_max[0]) return Q def baseline_empirical(model, demand_mu): """经验规则：库存<40%触发补货至80%""" M, T = model.M, model.T Q = np.zeros((M, T)) I_prev = model.I0 V = model.V for t in range(T): arrival = np.sum(Q[:, t]) * (1 - model.theta) I_after = I_prev + arrival - demand_mu if I_after < 0.4 * V: target = 0.8 * V order = min(target - I_after, model.Q_max[0]) Q[0, t] = max(0, order) I_prev = max(0, I_after) return Q def baseline_deterministic_opt(model, demand_mu): """确定性优化：以需求均值求解""" M, T = model.M, model.T Q = np.zeros((M, T)) I_prev = model.I0 target_inventory = demand_mu * 3 for t in range(T): arrival = np.sum(Q[:, t]) * (1 - model.theta) I_after = I_prev + arrival - demand_mu if I_after < target_inventory: gap = target_inventory - I_after order = min(gap, model.Q_max[0]) Q[0, t] = max(0, order) I_prev = max(0, I_after) return Q # ═══════════════════════════════════════════════════════════ # 3. GA-SA 混合启发式求解器 # ═══════════════════════════════════════════════════════════ class GASA_Solver: """遗传算法 + 模拟退火 混合启发式求解器""" def __init__(self, model, scenarios, probs=None, pop_size=150, max_gen=200, p_cross=0.9, p_mut=0.05, elite_rate=0.1, sa_iters=30, early_stop=20, verbose=True): self.model = model self.scenarios = scenarios self.probs = probs if probs is not None else np.ones(scenarios.shape[1])/scenarios.shape[1] self.M = model.M self.T = model.T self.pop_size = pop_size self.max_gen = max_gen self.p_cross = p_cross self.p_mut = p_mut self.elite_rate = elite_rate self.sa_iters = sa_iters self.early_stop = early_stop self.verbose = verbose # 自适应惩罚参数 self.penalty_lambda = 1.0 self.penalty_delta = 10.0 / max_gen def _init_population(self): """初始化种群：随机生成运输批量矩阵""" pop = [] for _ in range(self.pop_size): Q = np.zeros((self.M, self.T)) for i in range(self.M): # 每个补给源随机生成整数运输量 Q[i, :] = np.random.randint(0, int(self.model.Q_max[i]) + 1, self.T) pop.append(Q) return pop def _fitness(self, Q): """适应度函数 = 1/(总成本 + 惩罚)""" total, _, _, _, _ = self.model.compute_cost_scenario(Q, self.scenarios, self.probs) # 约束违反惩罚 penalty = 0 I_prev = self.model.I0 for t in range(self.T): arrival = np.sum(Q[:, t]) * (1 - self.model.theta) # 需求满足率惩罚（简化：检查每个场景的平均缺货） avg_shortage = np.mean([max(0, self.scenarios[t, s] - (I_prev + arrival)) for s in range(self.scenarios.shape[1])]) if avg_shortage > 0: penalty += self.model.Cs * avg_shortage * 10 I_after = max(0, I_prev + arrival - np.mean(self.scenarios[t, :])) # 库容惩罚 if I_after > self.model.V: penalty += 1e5 * (I_after - self.model.V) I_prev = I_after total_with_penalty = total + self.penalty_lambda * penalty return 1.0 / max(total_with_penalty, 1.0) def _tournament_select(self, pop, fitnesses, k=3): """锦标赛选择""" n = len(pop) selected = [] for _ in range(n): idx = np.random.choice(n, k, replace=False) winner = idx[np.argmax([fitnesses[i] for i in idx])] selected.append(deepcopy(pop[winner])) return selected def _crossover(self, p1, p2): """均匀交叉""" if np.random.random() > self.p_cross: return deepcopy(p1), deepcopy(p2) c1, c2 = deepcopy(p1), deepcopy(p2) mask = np.random.random((self.M, self.T)) < 0.5 c1[mask] = p2[mask] c2[mask] = p1[mask] return c1, c2 def _mutate(self, Q): """多项式变异""" Q_mut = deepcopy(Q) for i in range(self.M): for t in range(self.T): if np.random.random() < self.p_mut: delta = np.random.randint(-50, 51) Q_mut[i, t] = np.clip(Q_mut[i, t] + delta, 0, self.model.Q_max[i]) return Q_mut def _sa_local_search(self, Q): """模拟退火局部搜索""" current = deepcopy(Q) current_fit = self._fitness(current) best_Q = deepcopy(current) best_fit = current_fit # 初始温度设为种群目标函数标准差的10倍 T_k = 50000.0 for _ in range(self.sa_iters): # 邻域扰动：随机选择1-2个连续周期，重新分配运输量 neighbor = deepcopy(current) i = np.random.randint(self.M) t_a = np.random.randint(0, self.T - 1) t_b = min(t_a + np.random.randint(1, 4), self.T) total_ship = np.sum(current[i, t_a:t_b]) # 在窗口内随机重分配 if total_ship > 0: weights = np.random.random(t_b - t_a) weights /= weights.sum() neighbor[i, t_a:t_b] = np.clip( (total_ship * weights).astype(int), 0, self.model.Q_max[i] ) neighbor_fit = self._fitness(neighbor) delta = neighbor_fit - current_fit if delta > 0 or np.random.random() < np.exp(delta / max(T_k, 1e-6)): current = neighbor current_fit = neighbor_fit if current_fit > best_fit: best_Q = deepcopy(current) best_fit = current_fit T_k *= 0.95 return best_Q def solve(self): """运行GA-SA求解""" if self.verbose: print(f" GA-SA: pop={self.pop_size}, gen={self.max_gen}, SA_iters={self.sa_iters}") pop = self._init_population() best_Q = None best_fit = -np.inf no_improve = 0 for gen in range(self.max_gen): fitnesses = [self._fitness(ind) for ind in pop] # 更新最优 gen_best_idx = np.argmax(fitnesses) if fitnesses[gen_best_idx] > best_fit: best_fit = fitnesses[gen_best_idx] best_Q = deepcopy(pop[gen_best_idx]) no_improve = 0 else: no_improve += 1 # 精英保留 n_elite = max(1, int(self.pop_size * self.elite_rate)) elite_idx = np.argsort(fitnesses)[-n_elite:] elites = [deepcopy(pop[i]) for i in elite_idx] # 选择、交叉、变异 selected = self._tournament_select(pop, np.array(fitnesses)) offspring = [] for j in range(0, self.pop_size, 2): c1, c2 = self._crossover(selected[j], selected[min(j+1, self.pop_size-1)]) offspring.append(self._mutate(c1)) offspring.append(self._mutate(c2)) offspring = offspring[:self.pop_size] # 精英替换 for j in range(n_elite): offspring[j] = elites[j] # SA局部搜索精英个体 for j in range(n_elite): if np.random.random() < 0.5: offspring[j] = self._sa_local_search(offspring[j]) pop = offspring self.penalty_lambda += self.penalty_delta if self.verbose and gen % 50 == 0: total, h, s, t, l = self.model.compute_cost_scenario( best_Q, self.scenarios, self.probs ) sl = self.model.service_level(best_Q, self.scenarios, self.probs) print(f" Gen {gen:3d}: Z={total/1e4:.1f}万, SL={sl*100:.1f}%, λ={self.penalty_lambda:.1f}") if no_improve >= self.early_stop: if self.verbose: print(f" Early stop at gen {gen}") break total, h, s, t, l = self.model.compute_cost_scenario(best_Q, self.scenarios, self.probs) sl = self.model.service_level(best_Q, self.scenarios, self.probs) return best_Q, total, h, s, t, l, sl # ═══════════════════════════════════════════════════════════ # 4. MPC 滚动时域优化 # ═══════════════════════════════════════════════════════════ def mpc_rolling_horizon(model, scenarios, horizon=7, use_cluster=True): """ MPC滚动时域优化 每个周期仅优化未来H个周期，只执行第一个周期决策 """ M, T = model.M, model.T Q_full = np.zeros((M, T)) for t0 in range(T): h_end = min(t0 + horizon, T) H = h_end - t0 if H <= 0: break # 当前库存 I_current = model.I0 for tt in range(t0): I_current += np.sum(Q_full[:, tt]) * (1 - model.theta) I_current = max(0, I_current - np.mean(scenarios[tt, :])) # 为剩余H个周期优化（简化为贪心） sub_model = deepcopy(model) sub_model.T = H sub_model.I0 = I_current sub_scenarios = scenarios[t0:h_end, :] # 快速求解：使用贪心策略 Q_sub = np.zeros((M, H)) I_prev = I_current for h in range(H): # 贪心：补货至目标库存水平 target = np.mean(sub_scenarios[h, :]) * 3 gap = target - I_prev if gap > 0: for i in range(M): if gap <= 0: break order = min(gap, model.Q_max[i], gap * 0.7 if i == 0 else gap) Q_sub[i, h] = max(0, order) gap -= order I_prev = max(0, I_prev + np.sum(Q_sub[:, h]) * (1 - model.theta) - np.mean(sub_scenarios[h, :])) # 只执行第一个周期决策 Q_full[:, t0] = Q_sub[:, 0] return Q_full # ═══════════════════════════════════════════════════════════ # 5. 主程序 # ═══════════════════════════════════════════════════════════ def main(): print("=" * 70) print(" 成品油油库库存-运输协同优化求解器") print("=" * 70) # 加载数据 data_path = os.path.join(os.path.dirname(__file__), 'inventory_transport_data.xlsx') if not os.path.exists(data_path): data_path = '/home/ubuntu/chengzi-buluo/orange-family/docs/inventory_transport_data.xlsx' print(f"\n📂 数据文件: {data_path}") try: data = load_data_from_excel(data_path) print(f" M={data['M']}补给源, T={data['T']}周期") print(f" V={data['V']:.0f}吨, I0={data['I0']:.0f}吨, α={data['alpha']}") except Exception as e: print(f" ⚠️ 无法加载Excel，使用默认参数: {e}") # 使用默认参数 data = { 'M': 2, 'T': 30, 'V': 12000, 'I0': 3500, 'theta': 0.003, 'alpha': 0.95, 'Ch': 12, 'Cs': 850, 'Cl': 7200, 'sources': { 0: {'name': '炼厂A', 'Q_max': 800, 'L_km': 185, 'Ct': 0.45}, 1: {'name': '中转库B', 'Q_max': 600, 'L_km': 220, 'Ct': 0.52}, }, 'demand': np.random.RandomState(42).normal(385, 82, 30).clip(154, 770), 'scenarios': np.random.RandomState(77).normal(385, 82, (30, 100)).clip(115, None), } data['clustered'] = None data['cluster_probs'] = None model = InventoryTransportModel(data) demand = data['demand'] scenarios = data['scenarios'] # 使用聚合场景加速 if data.get('clustered') is not None: use_scenarios = data['clustered'] use_probs = data['cluster_probs'] print(f" 使用聚合场景: K={use_scenarios.shape[1]}") else: use_scenarios = scenarios[:, :10] # 只用前10个场景 use_probs = None print(f" 使用场景子集: S=10") demand_mu = demand.mean() demand_std = demand.std() # ── 运行所有方法 ── results = [] methods = [ ("(s,S)策略", lambda: baseline_ss_policy(model, demand_mu, demand_std)), ("独立库存优化(固定批量)", lambda: baseline_fixed_batch(model, 200, 3)), ("确定性优化", lambda: baseline_deterministic_opt(model, demand_mu)), ("经验规则(40%-80%)", lambda: baseline_empirical(model, demand_mu)), ] print("\n" + "=" * 70) print(" 求解中...") print("=" * 70) for name, fn in methods: t0 = time.time() Q = fn() elapsed = time.time() - t0 total, h, s, t, l = model.compute_cost_scenario(Q, use_scenarios, use_probs) sl = model.service_level(Q, use_scenarios, use_probs) results.append((name, Q, total, h, s, t, l, sl, elapsed)) print(f" {name:25s}: Z={total/1e4:8.2f}万, SL={sl*100:5.1f}%, T={elapsed:.2f}s") # GA-SA print(f"\n --- GA-SA 启发式求解 ---") ga_solver = GASA_Solver(model, use_scenarios, use_probs, pop_size=100, max_gen=150, sa_iters=20, early_stop=30, verbose=True) t0 = time.time() Q_gasa, total_gasa, h_gasa, s_gasa, t_gasa, l_gasa, sl_gasa = ga_solver.solve() elapsed_gasa = time.time() - t0 results.append(("GA-SA(本文)", Q_gasa, total_gasa, h_gasa, s_gasa, t_gasa, l_gasa, sl_gasa, elapsed_gasa)) # MPC print(f"\n --- MPC 滚动时域 ---") t0 = time.time() Q_mpc = mpc_rolling_horizon(model, use_scenarios, horizon=7) elapsed_mpc = time.time() - t0 total_mpc, h_mpc, s_mpc, t_mpc, l_mpc = model.compute_cost_scenario(Q_mpc, use_scenarios, use_probs) sl_mpc = model.service_level(Q_mpc, use_scenarios, use_probs) results.append(("MPC滚动时域", Q_mpc, total_mpc, h_mpc, s_mpc, t_mpc, l_mpc, sl_mpc, elapsed_mpc)) print(f" {'MPC滚动时域':25s}: Z={total_mpc/1e4:8.2f}万, SL={sl_mpc*100:5.1f}%, T={elapsed_mpc:.2f}s") # ── 输出结果汇总 ── print("\n" + "=" * 70) print(" 结果汇总") print("=" * 70) print(f"{'方法':<25s} {'Z(万元)':>10s} {'库存':>8s} {'缺货':>8s} {'运输':>8s} {'损耗':>8s} {'SL%':>6s} {'时间s':>6s}") print("-" * 82) best_idx = np.argmin([r[2] for r in results]) for i, (name, Q, total, h, s, t, l, sl, elapsed) in enumerate(results): marker = " ←最优" if i == best_idx else "" print(f"{name+marker:<25s} {total/1e4:10.2f} {h/1e4:8.2f} {s/1e4:8.2f} {t/1e4:8.2f} {l/1e4:8.2f} {sl*100:6.1f} {elapsed:6.2f}") # 相对改进 print("\n--- 相对改进（vs 经验规则） ---") base_total = results[3][2] # 经验规则 base_sl = results[3][7] for name, Q, total, h, s, t, l, sl, elapsed in results: cost_improve = (base_total - total) / base_total * 100 sl_improve = sl - base_sl print(f" {name:25s}: 成本 {cost_improve:+6.1f}%, 满足率 {sl_improve*100:+5.1f}pp") # ── 输出最优方案的运输计划 ── best_name, best_Q, _, _, _, _, _, _, _ = results[best_idx] print(f"\n--- 最优方案 ({best_name}) 运输计划 (前15周期) ---") print(f"{'周期':>4s}", end="") for i in range(model.M): print(f" {'Q' + str(i+1) + '(吨)':>10s}", end="") print() for t in range(min(15, model.T)): print(f"{t+1:4d}", end="") for i in range(model.M): print(f" {best_Q[i,t]:10.0f}", end="") print() # ── 保存结果到Excel ── save_results_to_excel(results, model, data_path) return results def save_results_to_excel(results, model, data_path): """保存结果到Excel""" from openpyxl import Workbook from openpyxl.styles import Font, PatternFill, Alignment, Border, Side out_path = data_path.replace('.xlsx', '_results.xlsx') wb = Workbook() hf = Font(bold=True, size=10, color="FFFFFF") hl = PatternFill(start_color="2F5496", end_color="2F5496", fill_type="solid") gf = PatternFill(start_color="C6EFCE", end_color="C6EFCE", fill_type="solid") tb = Border(left=Side('thin'),right=Side('thin'),top=Side('thin'),bottom=Side('thin')) s1 = wb.active; s1.title = "结果汇总" headers = ["方法","Z(万元)","库存成本","缺货成本","运输成本","损耗成本","满足率%","求解时间s"] for j, h in enumerate(headers, 1): c = s1.cell(row=1, column=j, value=h); c.font=hf; c.fill=hl; c.alignment=Alignment(horizontal='center'); c.border=tb best_idx = np.argmin([r[2] for r in results]) for i, (name, Q, total, h, s, t, l, sl, elapsed) in enumerate(results): vals = [name, round(total/1e4,2), round(h/1e4,2), round(s/1e4,2), round(t/1e4,2), round(l/1e4,2), round(sl*100,1), round(elapsed,2)] for j, v in enumerate(vals, 1): c = s1.cell(row=i+2, column=j, value=v); c.border=tb; c.alignment=Alignment(horizontal='center') if i == best_idx: c.fill = gf for c in range(1, 9): s1.column_dimensions[chr(64+c)].width = 14 # 最优方案详细计划 best_name, best_Q, _, _, _, _, _, _, _ = results[best_idx] s2 = wb.create_sheet("最优运输计划") s2.cell(row=1, column=1, value="周期") for i in range(model.M): s2.cell(row=1, column=i+2, value=f"Q{i+1}(吨)") for t in range(model.T): s2.cell(row=t+2, column=1, value=t+1) for i in range(model.M): s2.cell(row=t+2, column=i+2, value=int(best_Q[i,t])) for c in range(1, model.M+2): s2.column_dimensions[chr(64+c)].width = 14 wb.save(out_path) print(f"\n📊 结果已保存: {out_path}") if __name__ == '__main__': main()