线性规划灵敏度分析——Python实现

灵敏度分析（Sensitivity Analysis）是线性规划的一个重要部分，用于研究在模型参数发生变化时，最优解和目标函数值的变化情况。它能够识别和评估参数变动对解的影响，从而帮助决策者了解模型的稳定性及其对不同条件变化的反应。例如，通过灵敏度分析，决策者可以确定在什么范围内，目标函数系数、约束条件的右端常数或系数的变化不会改变最优解。这对于实际决策至关重要，因为它提供了关于系统弹性的关键洞察，允许决策者预见潜在风险和机会，并在资源供应、市场需求或成本变动时调整策略。灵敏度分析不仅可以用于优化当前计划，还可以用于制定应对未来不确定性的预案，从而实现更加稳健和合理的决策。这种分析方法在生产规划、投资决策、供应链管理和项目管理等领域广泛应用，帮助企业和组织在复杂多变的环境中保持竞争优势。

一、灵敏度分析理论

灵敏度分析主要涉及对以下三类参数的变化的研究：

目标函数系数的变化：研究目标函数的系数（如利润或成本系数）变化时对最优解和目标值的影响。
约束条件右端常数的变化：研究资源供应量或需求量等约束条件右端常数变化时对最优解的影响。
约束系数的变化：研究约束条件中的系数（如资源消耗系数）变化时对最优解的影响。

通过灵敏度分析，决策者可以回答以下问题：在什么范围内，目标函数系数的变化不会改变最优解的决策变量值？如果资源供应量发生变化，最优解的生产计划应如何调整？约束条件中的参数发生变化时，现有的最优解是否仍然有效？
这里以线性规划问题的标准型展开讨论：

$$\begin{aligned} & \max z=\boldsymbol{C} \boldsymbol{X} \\ & \text { s. t. }\left\{\begin{array}{c} \boldsymbol{A} \boldsymbol{X}=\boldsymbol{b} \\ \boldsymbol{X} \geq \mathbf{0} \end{array}\right. \end{aligned}$$

一般我们会添加松扡变量形成等式，所以初始单纯形表中的系数矩阵 $A$ 可以分为 $(B, N)$ ，其中 $B$为基变量对应的系数矩阵， $N$ 为非基变量对应的系数矩阵。同理，变量可分为基变量 $X_B$ ，非基变量 $X_N$ 。基变量对应的目标函数系数为 $C_B$ ，非基变量对应的目标函数系数为 $C_N$ 。
单纯形法本质上是在对增广矩阵进行多次初等行变换，使得最后满足可行和最优性检验，且基变量对应的系数矩阵为单位阵 $E$ 。根据线性代数的知识，我们知道，原来基变量对应的系数矩阵为 $B$ ，想要最后成为一个单位阵，只需左乘上 $B^{-1}$ 。因此，我们做了那么多次迭代，其实就相当于将增广矩阵左乘上一个 $B^{-1}$ 。
初始的基变量为松驰变量对应的矩阵 $N=E$ ，因此初始基变量在最优单纯形表中对应的系数矩阵即为 $B^{-1} N=B^{-1} E=B^{-1}$ 。右端资源向量 $b$ ，最终变为 $B^{-1} b$ 。下表给出了达到最优单纯形表时各个参数的特征。

$C_B$	$X_B$	$b$	$C_B$	$C_N$	$\theta$
			$X_B$	$X_N$
$C_B$	$X_B$	$B^{-1} b$	$B^{-1} B$	$B^{-1} N$
	$\sigma$		$C_B-C_B B^{-1} B$	$C_N-C_B B^{-1} N$	$z=C_B B^{-1} b$

有了这些特征后，我们可以清晰知道某个参数变化会带来哪些影响。如当右端资源向量 $b$ 变化时，只会影响最终基变量的取值 $B^{-1} b$ 。

1.1 目标函数系数 $c_j$ 的灵敏度分析

目标函数系数的灵敏度分析，主要考察在某个变量的系数变化时，最优解是否发生变化。具体而言，通过计算允许的增加量和减少量，可以确定目标函数系数在某个范围内变化时，不会改变当前的最优基解。通过分析最优单纯性表各个参数的特征，可知，目标函数系数 $c_j$ 的变化需要分情况讨论，因为 $c_j$ 为基变量和非基变量时，带来的影响变化不一样。

• $c_j$ 对应的变量 $x_j$ 为非基变量时
  此时 $c_j$ 的改变，只会影响对应决策变量 $x_j$ (非基) 的检验数。为保证基变量不变，该非基变量检验数需保持非正，即变化后的 $\sigma_j^{\prime}=c_j+\Delta c_j-c_B B^{-1} c_j \leq 0$ ，于是有:

$$\Delta c_j \leq-\sigma_j$$

当 $c_j$ 的变化值满足上述条件时，原来的线性规划模型最优解不变。若 $c_j$ 变动超过了上述范围，利用单纯形法继续求解。

• $c_j$ 对应的变量 $x_j$ 为基变量时
  此时 $c_j$ 进入 $C_B$ ，会影响所有决策变量的检验数，因此要保证所有非基变量检验数为非正。有:

$$\max \left(\sigma_j / \overline{a_{r j}} \mid \overline{a_{r j}}>0\right) \leq \Delta c_r \leq \min \left(\sigma_j / \overline{a_{r j}} \mid \overline{a_{r j}}<0\right)$$

式中， $\overline{a_{r j}}$ 为 $c_r$ 对应的基变量在最优单纯形表中的第 $r$ 行，第 $j$ 个元素。

1.2 约束条件右端常数 $b$ 的灵敏度分析

约束条件右端常数的灵敏度分析，关注的是资源供应量或需求量的变化对最优解的影响。右端资源列向量发生的变化的话，最优解中的基变量取值 $B^{-1} b$ 肯定会发生变化，因此我们要保证的是依然满足非负要求。有:

$$\max \left(-\overline{b_i} / \overline{a_{i r}} \mid \overline{a_{i r}}>0\right) \leq \Delta b_r \leq \min \left(-\overline{b_i} / \overline{a_{i r}} \mid \overline{a_{i r}}<0\right)$$

式中， $\overline{a_{i r}}$ 为 $B^{-1}$ 中的第 $i$ 行，第 $r$ 个元素。
当 $b_r$ 在上述范围内变动时，基变量和原来一样，但取值发生了变化，新基变量取值为 $B^{-1}\left(b_r+\right.$ $\left.\Delta b_r\right)$. 若变动超过了上述范围，不满足非负约束，需采用对偶单纯形法继续迭代求解。这时只能求出单个 $b$ 的变动范围，若变动了多个，一般直接求新解 $B^{-1} b^{\prime}$ 即可。

1.3 约束系数$a_{i j}$ 的灵敏度分析

约束系数的灵敏度分析，研究的是约束条件中系数的变化对最优解的影响。约束系数 $a_{i j}$ 的变化不仅可能会影响检验数，还可能会改变 $B^{-1}$。

• 约束系数 $a_{i j}$ 对应 $x_j$ 为非基变量
  此时只会改变该非基变量的检脸数，保证其非正的倩况下，用待定系数法进行变动范围的求解。
• 约束系数 $a_{i j}$ 对应 $x_j$ 为基变量
  此时 $B^{-1}$ 发生变化，假设变为 $\bar{B}^{-1}$ ，最终基变量取值 $B^{-1} b$ 则变为 $B^{-1} b$ ，所有非基变量检验数 $C_N-C_B B^{-1} N$ 变为 $C_N-$ $C_B \bar{B}^{-1} N$。此时会产生下面四种情况:
  • $\bar{B}^{-1} b \geq 0, C_N-C_B \bar{B}^{-1} N \leq 0$, 即此时既满足非负要求，又满足最优检验，已达到最优解。
  • $\bar{B}^{-1} b \geq 0$ 但 $C_N-C_B \bar{B}^{-1} N$ 存在正数，即此时虽然满足非负要求，但还有非惎变量检检数大于 0 ，用单靿形法继续选代求解。
  • $\bar{B}^{-1} b$ 存在负数。而 $C_N-C_B \bar{B}^{-1} N \leq 0$, 即此时最优检验满足了, 但非负性不满足, 因此利用对偶单纯形法继续迭代求解.
  • $\bar{B}^{-1} b$ 存在负数，且 $C_N-C_B \bar{B}^{-1} N$ 存在正数，即此时最优性检检和非负性要求均为达到。故需引进人工变量，用人工变量法或者从头算起。

二、灵敏度分析-python实现

一奶制品加工厂用牛奶生产A1,A2两种奶制品，1桶牛奶可以在甲车间用12小时加工成3公斤A1，或者在乙车间用8小时加工成4公斤A2。根据市场需求，生产的A1,A2全部能售出，且每公斤A1获利24元，每公斤A2获利16元。现在加工厂每天能得到50桶牛奶的供应，每天正式工人总的劳动时间480小时，并且甲车间每天至多能加工100公斤A1，乙车间的加工能力没有限制。试为该厂制订一个生产计划，使每天获利最大，并进一步讨论以下3个附加问题：
1） 若用35元可以买到1桶牛奶，应否作这项投资？
2） 若可以聘用临时工人以增加劳动时间，付给临时工人的工资最多是每小时几元？
3） 由于市场需求变化，每公斤A1的获利增加到30元，应否改变生产计划？

import pulp
import numpy as np

def solve_lp():
    # 定义问题
    problem = pulp.LpProblem("Maximize_Profit", pulp.LpMaximize)

    # 定义决策变量
    A1 = pulp.LpVariable("A1", lowBound=0)
    A2 = pulp.LpVariable("A2", lowBound=0)

    # 目标函数
    problem += 72 * A1 + 64 * A2, "Total_Profit"

    # 约束条件
    problem += A1 + A2 <= 50, "Milk_Constraint"
    problem += 12 * A1 + 8 * A2 <= 480, "Labor_Constraint"
    problem += A1 <= 100, "A1_Constraint"

    # 求解问题
    problem.solve()

    # 输出结果
    print(f"Status: {pulp.LpStatus[problem.status]}")
    print(f"A1 Production: {A1.varValue:.2f} kg")
    print(f"A2 Production: {A2.varValue:.2f} kg")
    print(f"Total Profit: {pulp.value(problem.objective):.2f}")

    # 获取基础变量和约束
    basis_matrix = np.array([[1, 1], [12, 8]])
    basis_inverse = np.linalg.inv(basis_matrix)

    print("Basis Matrix:")
    print(basis_matrix)
    print("Inverse of Basis Matrix:")
    print(basis_inverse)

    return problem, basis_inverse

def sensitivity_analysis(problem, basis_inverse):
    # 提取对偶变量（影子价格）和松弛量
    milk_shadow_price = problem.constraints['Milk_Constraint'].pi
    labor_shadow_price = problem.constraints['Labor_Constraint'].pi

    # 当前约束右侧常数
    milk_rhs = 50
    labor_rhs = 480

    # 获取最优解
    A1_opt = pulp.value(problem.variablesDict()['A1'])
    A2_opt = pulp.value(problem.variablesDict()['A2'])
    A_opt = np.array([A1_opt, A2_opt])
     
    # 计算变化区间
    milk_change_lower = milk_rhs + max([-A_opt[j]/ basis_inverse[j,0] for j in range(basis_inverse.shape[1]) if basis_inverse[0, j] > 0])
    milk_change_upper = milk_rhs + min([-A_opt[j]/ basis_inverse[j,0] for j in range(basis_inverse.shape[1]) if basis_inverse[0, j] < 0])

    labor_change_lower = labor_rhs + max([-A_opt[j] / basis_inverse[j, 1] for j in range(basis_inverse.shape[1]) if basis_inverse[1, j] > 0])
    labor_change_upper = labor_rhs + min([-A_opt[j] / basis_inverse[j, 1] for j in range(basis_inverse.shape[1]) if basis_inverse[1, j] < 0])

    print("\nSensitivity Analysis:")
    print(f"Milk Constraint: Shadow Price = {milk_shadow_price:.2f}")
    print(f"Milk Supply Change Interval: [{milk_change_lower:.2f}, {milk_change_upper:.2f}]")
    print(f"Labor Constraint: Shadow Price = {labor_shadow_price:.2f}")
    print(f"Labor Supply Change Interval: [{labor_change_lower:.2f}, {labor_change_upper:.2f}]")

if __name__ == "__main__":
    problem, basis_inverse = solve_lp()
    sensitivity_analysis(problem, basis_inverse)

A1 Production: 20.00 kg
A2 Production: 30.00 kg
Total Profit: 3360.00

Sensitivity Analysis:
Milk Constraint: Shadow Price = 48.00
Milk Supply Change Interval: [40.00, 60.00]
Labor Constraint: Shadow Price = 2.00
Labor Supply Change Interval: [400.00, 600.00]

这个线性规划的最优解为A1=20，A2=30，最优值为z=3360，即用20桶牛奶生产A1, 30桶牛奶生产A2，可获最大利润3360元。输出中除了告诉我们问题的最优解和最优值以外，还有许多对分析结果有用的信息。 输出中Shadow Price给出这2种资源在最优解下“资源”增加1个单位时“效益”的增量：原料增加1个单位（1桶牛奶）时利润增长48（元），劳动时间增加1个单位（1小时）时利润增长2（元）。用影子价格的概念很容易回答附加问题(1）：用35元可以买到1桶牛奶，低于1桶牛奶的影子价格48，当然应该作这项投资。回答附加问题(2）：聘用临时工人以增加劳动时间，付给的工资低于劳动时间的影子价格才可以增加利润，所以工资最多是每小时2元。对每公斤A1的获利增加到30元，可用A1的价值系数进行灵敏度分析，也可以将A1的价值系数调整到30，重新求出其最优解，可知不改变最优生产计划。

三、案例分析

已知某生产计划的线性规划模型如下:

$$\text { s.t. }\left\{\begin{array}{c} \max z=x_1+5 x_2+3 x_3+4 x_4 \\ 5 x_1+3 x_2+x_3+2 x_4 \leq 800 \\ 3 x_1+4 x_2+3 x_3+4 x_4 \leq 1200 \\ x_1+4 x_2+5 x_3+3 x_4 \leq 1000 \\ x_j \geq 0, j=1,2,3,4 \end{array}\right.$$

$x_1, x_2, x_3, x_4$ 分别代表产品的产量, $x_5, x_6, x_7$ 分别是对应资源 $1,2,3$ 的松驰变量, 利用单纯形法得到的最优单纯形表如下:

$c_j$			1	5	3	4	0	0	0
$C_B$	$\boldsymbol{X}_B$	$b$	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$x_7$
0	$x_5$	100	0.25	0	-3.25	0	$\mathbf{1}$	0.25	-1
4	$x_4$	200	2	0	-2.00	$\mathbf{1}$	0	1	-1
5	$x_2$	100	-0.75	$\mathbf{1}$	2.75	0	0	-0.75	1
$z_j$			4.25	5	5.75	4	0	0.25	1
$c_j-z_j$			-3.25	0	-2.75	0	0	-0.25	-1

(1) 对 $c_4$ 进行灵敏度分析。
(2) 现在准备生产新产品, 已知该新产品耗用资源 1 的数量为 5 , 耗用资源 2 的数量为 4 , 耗用资源 3 的数量为 3 , 该新产品的单位利润为 9 , 那么生产该新产品是否有利? 为什么? 如果生产此产品有利, 新的最优生产方案是什么?

3.1 $c_4$ 灵敏度分析

$c_4$ 进行灵敏度分析见下面程序计算。

import numpy as np
from scipy.optimize import linprog

# Coefficients of the objective function (for maximization, we negate the coefficients)
c = [-1, -5, -3, -4, 0, 0, 0]

# Coefficients of the equality constraints
A_eq = np.array([
    [5, 3, 1, 2, 1, 0, 0],
    [3, 4, 3, 4, 0, 1, 0],
    [1, 4, 5, 3, 0, 0, 1]
])

# Right-hand side of the equality constraints
b_eq = [800, 1200, 1000]

# Bounds for the variables
x_bounds = (0, None)  # x_j >= 0 for all j

# Solve the linear programming problem using the HiGHS solver
res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=[x_bounds]*7, method='highs')

# Output the results
if res.success:
    print(f"Optimal value of objective function: {-res.fun}")  # Negate back to original maximization problem
    for i, x in enumerate(res.x[:4]):
        print(f"x{i+1} = {x}")
    for i, s in enumerate(res.x[4:], start=1):
        print(f"s{i} = {s}")

    # Identify the basic variables by finding non-zero values in the solution
    basic_indices = np.where(res.x > 1e-5)[0]  # Use a small threshold to account for numerical precision
    print("\nBasic Variable Indices:")
    print(basic_indices)

    # Extract the coefficients of the basic variables from c
    C_B = np.array([c[i] for i in basic_indices])
    print("\nCoefficients of Basic Variables (C_B):")
    print(C_B)

    # Form the basis matrix B using columns of A_eq corresponding to basic variables
    B = A_eq[:, basic_indices]
    print("\nBasis Matrix B:")
    print(B)

    # Compute the inverse of the basis matrix
    if B.shape[0] == B.shape[1]:  # Ensure it's square before inversion
        B_inv = np.linalg.inv(B)
        print("\nInverse of the Basis Matrix:")
        print(B_inv)
    else:
        print("\nThe basis matrix is not square, and its inverse cannot be computed directly.")

    # Compute the reduced costs (σ_j')
    reduced_costs = np.array(c) - C_B @ B_inv @ A_eq
    print("\nReduced Costs (σ'):")
    for i, rc in enumerate(reduced_costs):
        print(f"σ_{i+1}' = {rc}")

    # Compute the interval for Δc_4
    basis_inverse_row_2 = B_inv[1] @ A_eq
    valid_lower = [-reduced_costs[j] / basis_inverse_row_2[j] for j in range(len(basis_inverse_row_2)) if basis_inverse_row_2[j] > 0 and reduced_costs[j] != 0]
    valid_upper = [-reduced_costs[j] / basis_inverse_row_2[j] for j in range(len(basis_inverse_row_2)) if basis_inverse_row_2[j] < 0 and reduced_costs[j] != 0]
           
    c4_lower = max(valid_lower) - c[3] if valid_lower else float('-inf')
    c4_upper = min(valid_upper) - c[3] if valid_upper else float('inf')
    print(f"\nInterval for Δc_4: [{max(valid_lower)}, {min(valid_upper)}]")
    print(f"\nInterval for c_4: [{c4_lower}, {c4_upper}]")

else:
    print("No optimal solution found.")

Interval for Δc_4: [-0.25, 1.0]
Interval for c_4: [3.75, 5.0]

3.2 新产品是否生产

添加一个新列重新计算其最优解，可回答问题（2）生产新产品。

import numpy as np
from scipy.optimize import linprog

# Coefficients of the objective function (for maximization, we negate the coefficients)
c = [-1, -5, -3, -4, 0, 0, 0, -9]  # Including profit of the new product

# Coefficients of the equality constraints
A_eq = np.array([
    [5, 3, 1, 2, 1, 0, 0, 5],     
    [3, 4, 3, 4, 0, 1, 0, 4],    
    [1, 4, 5, 3, 0, 0, 1, 3],    
])

# Right-hand side of the equality constraints
b_eq = [800, 1200, 1000]  # The last entry is 0 for resource availability for the new product

# Bounds for the variables
x_bounds = (0, None)  # x_j >= 0 for all j

# Solve the linear programming problem using the HiGHS solver
res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=[x_bounds]*8, method='highs')

# Output the results
if res.success:
    # Print optimal solution
    print(f"Optimal value of objective function: {-res.fun}")  # Negate back to original maximization problem
    for i, x in enumerate(res.x[:4]):
        print(f"x{i+1} = {x}")
    for i, s in enumerate(res.x[4:7], start=1):
        print(f"s{i} = {s}")
    print(f"Production quantity of new product: {res.x[7]}")

    # Identify the basic variables by finding non-zero values in the solution
    basic_indices = np.where(res.x > 1e-5)[0]  # Use a small threshold to account for numerical precision
    print("\nBasic Variable Indices:")
    print(basic_indices)

    # Extract the coefficients of the basic variables from c
    C_B = np.array([c[i] for i in basic_indices])
    print("\nCoefficients of Basic Variables (C_B):")
    print(C_B)

    # Form the basis matrix B using columns of A_eq corresponding to basic variables
    B = A_eq[:, basic_indices]
    print("\nBasis Matrix B:")
    print(B)

    # Compute the inverse of the basis matrix
    if B.shape[0] == B.shape[1]:  # Ensure it's square before inversion
        B_inv = np.linalg.inv(B)
        print("\nInverse of the Basis Matrix:")
        print(B_inv)
    else:
        print("\nThe basis matrix is not square, and its inverse cannot be computed directly.")

    # Compute the reduced costs (σ_j')
    reduced_costs = np.array(c) - C_B @ B_inv @ A_eq
    print("\nReduced Costs (σ'):")
    for i, rc in enumerate(reduced_costs):
        print(f"σ_{i+1}' = {rc}")

    # Determine if producing the new product is profitable
    if res.x[7] > 0:
        print("\nProducing the new product is profitable.")
    else:
        print("\nProducing the new product is not profitable.")

else:
    print("No optimal solution found.")

总结

灵敏度分析是线性规划中一个重要的工具，帮助决策者理解模型的稳定性及其对参数变化的反应。它能够识别和评估参数变动对解的影响，从而使决策者能够在不确定性环境中做出更加稳健的决策。通过灵敏度分析，决策者可以了解在什么范围内目标函数系数、约束条件的右端常数或系数的变化不会改变最优解。这对于资源配置、成本控制和风险管理至关重要。在实际应用中，灵敏度分析能够帮助企业优化资源配置，调整生产计划，提高整体效益。例如，企业可以通过分析原材料价格或市场需求的变化，预测其对生产计划和利润的影响，从而制定应对策略。在现代复杂的商业环境中，灵敏度分析的应用前景广阔，具有重要的理论和实践价值。它不仅可以帮助企业在面对市场波动、政策变化和技术革新等不确定因素时保持竞争力，还可以为政府部门在公共政策制定、预算分配和社会经济规划等方面提供科学依据。

参考文献

1. 一个实例理解Lingo的灵敏性分析
2. 【管理运筹学】第 4 章——线性规划问题的灵敏度分析
3. 线性规划问题的灵敏度分析