第 5 章：标准形式与灵敏度分析 —— 统一度量衡与参数敏感度¶

# 示例：max → min
original_max = "max z = 5x1 + 3x2"
converted_min = "min z' = -5x1 - 3x2  (其中 z' = -z)"

print(f"原始: {original_max}")
print(f"转换: {converted_min}")
print(f"最优时: max z = -{converted_min.split('=')[1].strip()}")

def convert_to_standard():
    """
    将以下问题化为标准形式：
    max z = 5x1 + 3x2
    s.t. 2x1 + x2 <= 10
         x1 + 2x2 >= 8
         x1, x2 >= 0
    """
    print("原始问题:")
    print("  max z = 5x1 + 3x2")
    print("  s.t. 2x1 + x2 <= 10")
    print("        x1 + 2x2 >= 8")
    print("        x1, x2 >= 0")

    print("\n转换步骤:")

    # Step 1: max → min
    print("  1. max → min: 令 z' = -z")
    print("     min z' = -5x1 - 3x2")

    # Step 2: 松弛变量
    print("  2. 第一个约束 + 松弛变量 s1:")
    print("     2x1 + x2 + s1 = 10,  s1 >= 0")

    # Step 3: 剩余变量
    print("  3. 第二个约束 - 剩余变量 s2:")
    print("     x1 + 2x2 - s2 = 8,  s2 >= 0")

    print("\n标准形式:")
    print("  min z' = -5x1 - 3x2 + 0*s1 + 0*s2")
    print("  s.t. 2x1 + x2 + s1       = 10")
    print("        x1 + 2x2 - s2       = 8")
    print("        x1, x2, s1, s2 >= 0")

convert_to_standard()

def basic_solutions_demo():
    """
    标准形式的基与基本解

    min z' = -5x1 - 3x2 + 0*s1 + 0*s2
    s.t. 2x1 + x2 + s1       = 10
          x1 + 2x2 - s2       = 8
          x1, x2, s1, s2 >= 0

    变量顺序: [x1, x2, s1, s2]
    """
    print("约束矩阵 A =")
    A = np.array([[2, 1, 1, 0], [1, 2, 0, -1]])
    print(A)

    print("\n可能的基组合（选2个作为基变量）:")

    bases = [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]

    for base in bases:
        non_base = [i for i in range(4) if i not in base]
        print(f"\n基变量: {[chr(120+int(b)) for b in base]} (索引 {base})")
        print(f"非基变量: {[chr(120+int(nb)) for nb in non_base]} = 0")

        # 提取基矩阵
        B = A[:, base]
        b = np.array([10, 8])

        # 判断基矩阵是否可逆
        det = np.linalg.det(B)
        if abs(det) < 1e-10:
            print(f"  基矩阵奇异，跳过")
            continue

        # 求解
        x_base = np.linalg.solve(B, b)

        if np.all(x_base >= -1e-10):
            z = -5 * x_base[0] - 3 * x_base[1]
            print(f"  解: {[chr(120+int(b)) + '=' + str(round(x_base[i], 2)) for i, b in enumerate(base)]}")
            print(f"      z' = {z:.2f}, z = {-z:.2f} ✓ 可行")
        else:
            print(f"  解包含负值，不可实施")

basic_solutions_demo()

def sensitivity_c1(c1, c2=3):
    """
    分析 c1 变化时，最优解是否保持不变

    最优基: [x1, x2] = (4, 2)
    最优基矩阵 B = [[2, 1], [1, 2]]
    """
    B = np.array([[2, 1], [1, 2]])
    B_inv = np.linalg.inv(B)

    # 最优性条件: 检验数 sigma_N <= 0 (对于 min 问题)
    # 非基变量 s1 的检验数: c_B * B^(-1) * a_s1 - c_s1
    # 非基变量 s2 的检验数: c_B * B^(-1) * a_s2 - c_s2

    c_B = np.array([c1, c2])  # 基变量成本

    # a_s1 = [1, 0], c_s1 = 0
    sigma_s1 = np.dot(c_B, B_inv[:, 0]) - 0

    # a_s2 = [0, -1], c_s2 = 0
    sigma_s2 = np.dot(c_B, B_inv[:, 1]) - 0

    return sigma_s1, sigma_s2

print("灵敏度分析：c1 变化对最优性的影响")
print("=" * 55)
print(f"{'c1范围':<20} {'σ_s1':<12} {'σ_s2':<12} {'最优基稳定?'}")
print("-" * 55)

for c1 in [2, 3, 4, 5, 6, 7, 8]:
    sigma_s1, sigma_s2 = sensitivity_c1(c1)
    stable = "是" if sigma_s1 <= 0 and sigma_s2 <= 0 else "否"
    print(f"{c1:<20} {sigma_s1:<12.4f} {sigma_s2:<12.4f} {stable}")

def shadow_price_analysis():
    """
    影子价格：右端项每增加1单位，最优目标函数值的增量

    B^(-1) 的第一列 = [2/3, -1/3]
    """
    B = np.array([[2, 1], [1, 2]])
    B_inv = np.linalg.inv(B)

    # 最优解时，影子价格 y* = c_B * B^(-1)
    c_B = np.array([5, 3])
    y_star = np.dot(c_B, B_inv)

    print(f"基矩阵 B 的逆:\n{B_inv}")
    print(f"\n影子价格 y* = c_B * B^(-1) = {y_star}")
    print(f"\n解读:")
    print(f"  - 增加 1 单位牛奶（右端项1）: 利润增加 {y_star[0]:.4f} 元")
    print(f"  - 增加 1 单位茶叶（右端项2）: 利润增加 {y_star[1]:.4f} 元")
    print(f"\n这意味着：牛奶比茶叶更值钱（每单位资源带来的边际利润）")

shadow_price_analysis()