《交通运输系统工程》课程辅助教学平台

import numpy as np

def uncertainty_decision_analysis(payoff_matrix, beta=0.6):
    """
    不确定型决策五大准则计算引擎
    输入: payoff_matrix (各方案在不同状态下的益损值矩阵), beta (折中系数)
    输出: 5种准则下的最优方案索引
    """
    matrix = np.array(payoff_matrix)
    
    # 1. 悲观准则 (Max-Min, 小中取大)
    pes_min_vals = np.min(matrix, axis=1)
    best_pes = np.argmax(pes_min_vals)
    
    # 2. 乐观准则 (Max-Max, 大中取大)
    opt_max_vals = np.max(matrix, axis=1)
    best_opt = np.argmax(opt_max_vals)
    
    # 3. 折中准则 (Hurwicz准则)
    hurwicz_vals = beta * opt_max_vals + (1 - beta) * pes_min_vals
    best_hur = np.argmax(hurwicz_vals)
    
    # 4. 等可能准则 (Laplace准则)
    laplace_vals = np.mean(matrix, axis=1)
    best_lap = np.argmax(laplace_vals)
    
    # 5. 遗憾准则 (Minimax Regret)
    # 首先求每一列(每种状态)的最大益损值
    col_max = np.max(matrix, axis=0) 
    # 构造后悔值矩阵 (某状态最大收益 - 当前方案收益)
    regret_matrix = col_max - matrix
    # 求每个方案的最大后悔值
    max_regrets = np.max(regret_matrix, axis=1)
    # 在最大后悔值中取最小
    best_reg = np.argmin(max_regrets)
    
    return best_pes, best_opt, best_hur, best_lap, best_reg

# ================= 运行例题 9-1 数据 =================
payoff_data = [
    [110, 50, 40], # A1: 大型仓库
    [70,  80, 50], # A2: 中型仓库
    [30,  60, 70]  # A3: 小型仓库
]
res = uncertainty_decision_analysis(payoff_data, beta=0.6)
alts = ["A1(大型)", "A2(中型)", "A3(小型)"]
print("悲观准则最优:", alts[res[0]])
print("乐观准则最优:", alts[res[1]])
print("折中准则最优:", alts[res[2]])
print("等可能最优  :", alts[res[3]])
print("遗憾准则最优:", alts[res[4]])

# 单级决策树 Python 求解示例
def solve_single_level(p_t1, p_t2, p_t3, p_t4):
    # 各方案在不同状态下的益损值表
    payoffs = {
        'A1(大型)': [180, 90, -30, -60],
        'A2(中型)': [120, 60, 20, -10],
        'A3(租赁)': [80, 30, 40, 10]
    }
    probs = [p_t1, p_t2, p_t3, p_t4]
    
    expected_values = {}
    for plan, values in payoffs.items():
        # 计算期望值: sum(收益 * 概率)
        ev = sum(v * p for v, p in zip(values, probs))
        expected_values[plan] = ev
        
    # 取最大值对应的方案
    best_plan = max(expected_values, key=expected_values.get)
    return best_plan, expected_values

# 运行教材默认概率
best, evs = solve_single_level(0.4, 0.2, 0.3, 0.1)
print(f"各方案期望值: {evs}")
print(f"最优方案: {best}")

# 多阶段决策树 Python 逆向推导示例
def solve_multi_level(p_rain):
    p_no_rain = 1 - p_rain
    
    # ---------------- 阶段二 (节点 Z) ----------------
    # Z1 (紧急加班) 的预期损失
    ev_z1 = 0.5 * (-30000) + 0.3 * (-22500) + 0.2 * (-15000)
    # Z2 (正常施工) 的固定损失
    ev_z2 = -25000
    
    # Z节点决策: 取损失最小(数值最大)
    best_z_val = max(ev_z1, ev_z2) 
    
    # ---------------- 阶段一 (节点 Y 和 X) ----------------
    # Y节点状态期望值
    ev_y = p_no_rain * 0 + p_rain * best_z_val
    
    # X节点方案
    ev_x1 = -15000 #直接加班
    ev_x2 = ev_y   # 正常施工
    
    # X节点决策
    if ev_x1 > ev_x2:
        return "选择方案 X1(直接加班)", ev_x1, ev_x2
    else:
        return "选择方案 X2(正常施工)", ev_x1, ev_x2

decision, x1_cost, x2_cost = solve_multi_level(0.7)
print(f"X1期望损失: {x1_cost}, X2期望损失: {x2_cost}")
print(f"最终决策: {decision}")

import numpy as np

def gm11_predict(x0, forecast_steps=1):
    """灰色预测模型 GM(1,1)"""
    n = len(x0)
    # 1. 累加生成 (1-AGO)
    x1 = np.cumsum(x0) 
    # 2. 构造背景值及 B 矩阵、Y 向量
    z1 = 0.5 * (x1[:-1] + x1[1:]) 
    B = np.vstack((-z1, np.ones_like(z1))).T
    Y = x0[1:].reshape(-1, 1)
    # 3. 最小二乘法求参数 a, u
    a, u = np.linalg.inv(B.T @ B) @ B.T @ Y 
    a, u = a[0], u[0]
    # 4. 时间响应函数
    def response(k): 
        return (x0[0] - u/a) * np.exp(-a * k) + u/a 
    # 5. 累减还原
    pred_x1 = np.array([response(i) for i in range(n + forecast_steps)]) 
    pred_x0 = np.concatenate(([x0[0]], np.diff(pred_x1)))
    return a, u, pred_x0

# ================= 运行例题 6-8 数据 =================
traffic_data = np.array([6986, 8013, 9353, 10821, 12304, 13755])
a, u, results = gm11_predict(traffic_data, forecast_steps=1)
print(f"发展系数 a: {a:.5f}\n灰色作用量 u: {u:.3f}\n2025年预测值: {results[-1]:.0f} pcu/d")

import numpy as np

def ahp_analyze(matrix):
    """层次分析法核心计算：特征根法"""
    n = matrix.shape[0]
    eigenvalues, eigenvectors = np.linalg.eig(matrix)
    max_idx = np.argmax(np.real(eigenvalues))
    lambda_max = np.real(eigenvalues[max_idx])
    w = np.real(eigenvectors[:, max_idx])
    weights = w / np.sum(w)
    
    CI = (lambda_max - n) / (n - 1)
    RI_dict = {1:0, 2:0, 3:0.52, 4:0.89, 5:1.12, 6:1.26, 7:1.36, 8:1.41, 9:1.46, 10:1.49}
    RI = RI_dict.get(n, 1.54)
    CR = (CI / RI) if (n > 2 and RI > 0) else 0.0
    return weights, lambda_max, CI, CR, CR < 0.1

# ================= 运行例题 8-9 的 C-G 矩阵 =================
matrix = [[1, 3, 7, 9], [1/3, 1, 5, 7], [1/7, 1/5, 1, 3], [1/9, 1/7, 1/3, 1]]
w, l_max, ci, cr, valid = ahp_analyze(np.array(matrix))
print(f"权重向量 W: {np.round(w, 4)}\n一致性比例 CR: {cr:.4f}")

import numpy as np

def group_ahp_synthesis(expert_weights, individual_vectors):
    lambdas = np.array(expert_weights)
    lambdas = lambdas / np.sum(lambdas) # 归一化
    final_weights = np.dot(lambdas, individual_vectors)
    return final_weights

# 专家权重系数
lambdas = [0.15, 0.15, 0.1, 0.1, 0.05, 0.1, 0.1, 0.05, 0.15, 0.05]
expert_1_w = np.array([0.0833, 0.0417, 0.1125, 0.0375, 0.1125, 0.1125, 
                       0.0536, 0.1071, 0.1071, 0.1071, 0.0938, 0.0313])
# 综合 10 位专家的向量计算 final_w
final_w = group_ahp_synthesis(lambdas, all_expert_vectors)
print("最终权重:", np.round(final_w, 4))

import numpy as np

def fuzzy_comprehensive_evaluation(weights, relation_matrix, scales):
    """
    模糊综合评价核心计算引擎
    输入: 
      weights: 一维权重向量 A (要素权重)
      relation_matrix: 隶属度矩阵 R (方案的模糊评价分布)
      scales: 评价等级打分尺度 V
    输出:
      B: 模糊变换后的综合评价向量
      d: 标量化的可行度得分
    """
    A = np.array(weights)
    R = np.array(relation_matrix)
    V = np.array(scales)
    
    # 1. 模糊变换运算 B = A * R (在此使用普通矩阵乘法计算模糊算子)
    B = np.dot(A, R)
    
    # 2. 计算方案可行度得分 d = B * V^T
    d = np.dot(B, V.T)
    
    return B, d

# ================= 运行例题 8-11 数据 =================
A_weights = [0.2, 0.3, 0.5]
V_scales = [5, 3, 1]

# 方案的隶属度矩阵
R_dict = {
    '方案A': [[0.7, 0.2, 0.1], [0.1, 0.2, 0.7], [0.3, 0.6, 0.1]],
    '方案B': [[0.3, 0.6, 0.1], [1.0, 0.0, 0.0], [0.7, 0.3, 0.0]],
    '方案C': [[0.1, 0.4, 0.5], [1.0, 0.0, 0.0], [0.1, 0.3, 0.6]]
}

scores = {}
for name, R in R_dict.items():
    B_res, d_res = fuzzy_comprehensive_evaluation(A_weights, R, V_scales)
    scores[name] = d_res
    print(f"{name} -> 模糊向量 B: {np.round(B_res, 2)}, 可行度得分 d: {d_res:.2f}")

best_scheme = max(scores, key=scores.get)
print(f"\n最佳选择为: {best_scheme}")

import numpy as np

def calculate_evpi(cost_matrix, probabilities):
    """
    计算完全信息的价值 (针对成本最小化问题)
    """
    costs = np.array(cost_matrix)
    probs = np.array(probabilities)
    
    # 1. 计算无信息期望成本 EMV
    emv_array = np.dot(costs, probs)
    min_emv = np.min(emv_array)
    
    # 2. 计算具备完全信息时的期望成本 EPPI
    min_costs_per_state = np.min(costs, axis=0)
    eppi = np.dot(min_costs_per_state, probs)
    
    # 3. 完全信息的价值 EVPI
    evpi = min_emv - eppi
    
    return min_emv, eppi, evpi

# 运行例 9-5 数据
costs = [[500, 800], [700, 600]]
probs = [0.6, 0.4]

emv, eppi, evpi = calculate_evpi(costs, probs)
print(f"完全信息价值 EVPI: {evpi}")

# 例 9-6 抽样信息的价值计算
# 1. 不测试直接投放的期望收益
emv_direct = 0.75 * 1000 + 0.25 * (-180)

# 2. 测试分支逆向归纳计算
# 若测试结果为好 (概率 0.85)
emv_test_good_invest = 0.9 * 1000 + 0.1 * (-180) # = 882
emv_test_good_no_invest = 0
# 剪枝：选大的
best_if_good = max(emv_test_good_invest, emv_test_good_no_invest) # 882

# 若测试结果为差 (概率 0.15)
emv_test_bad_invest = 0.15 * 1000 + 0.85 * (-180) # = -3
emv_test_bad_no_invest = 0
# 剪枝：选大的
best_if_bad = max(emv_test_bad_invest, emv_test_bad_no_invest) # 0

# 3. 综合测试节点的期望值
emv_test = 0.85 * best_if_good + 0.15 * best_if_bad # = 749.7

# 4. 计算 EVSI
evsi = emv_test - emv_direct

print(f"不测试最大期望: {emv_direct} 万")
print(f"测试分支总期望: {emv_test} 万")
print(f"抽样信息价值 EVSI: {evsi:.1f} 万")

# 例 9-7 效用理论与决策

def utility(x, alpha, x_min=-1000, x_max=1400):
    """根据风险偏好系数计算归一化效用值"""
    if x <= x_min: return 0.0
    if x >= x_max: return 1.0
    return ((x - x_min) / (x_max - x_min)) ** (1.0 / alpha)

# 方案1：高速公路
v_high_good = 200 * 10 - 600  # = 1400
v_high_bad = -40 * 10 - 600   # = -1000
# 方案2：二级公路
v_low_good = 80 * 10 - 280    # = 520
v_low_bad = 60 * 10 - 280     # = 320

p_good = 0.7
p_bad = 0.3

# 测试：如果是保守型决策者 (alpha = 2.0)
alpha = 2.0

eu_high = p_good * utility(v_high_good, alpha) + p_bad * utility(v_high_bad, alpha)
eu_low = p_good * utility(v_low_good, alpha) + p_bad * utility(v_low_bad, alpha)

print(f"高速公路期望效用 EU: {eu_high:.3f}")
print(f"二级公路期望效用 EU: {eu_low:.3f}")

if eu_high > eu_low:
    print("决策：建设高速公路")
else:
    print("决策：建设二级公路 (保守选择)")