线性规划灵敏度分析精解

灵敏度分析(Sensitivity Analysis),又称后优化分析,是在获得线性规划问题的最优解之后,对相关参数变化进行系统研究的过程,以确定解的稳定性和鲁棒性。通过灵敏度分析,我们可以评估不同参数变化对最优解和目标值的影响,从而了解模型的响应性和可靠性。具体而言,这些参数包括目标函数的系数(例如利润或成本系数)、约束条件的右端项(如资源的供应量或需求量)以及约束系数矩阵的元素(即资源的消耗系数或生产技术系数)。灵敏度分析的核心目的是探讨这些参数在发生变化时,最优解是否依然保持不变或者会发生何种变化,进而帮助决策者更好地应对不确定性,优化资源配置,并在面对外部条件变化时,及时调整策略,以确保决策的有效性和最优性。

灵敏度分析思路 灵敏度分析内容

一、灵敏度分析的内容

讨论线性规划问题时,尝尝假定A, b, C都是常数。但实际上这些系数往往是估计值和预测值。如市场条件一变,C值就会变化;A往往是因工艺技术条件的改变而改变;b是根据资源投入后的经济效果决定的一种决策选择。因此提出这样两个问题:

当这些系数有一个或几个发生变化时,已求得的线性规划问题的最优解会有什么变化;
或者这些系数在什么范围内变化时,线性规划问题的最优解或最优基不变。

f(x)=cBB1b(cBB1NcN)xN

灵敏度分析 xB xN
xB Im B1N B1b
f 0 cNcBB1N cBB1b

‌灵敏度分析通过揭示参数变化对决策结果的潜在影响,使得决策者能够在变化的环境中保持灵活性和适应性,从而做出更加明智和有效的决策。

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

目标函数系数表示决策变量对目标函数的贡献程度,通过改变目标函数系数可以分析目标函数值的变动情况。当目标函数系数发生较大变动时,可能导致最优解的决策变量发生改变。
企业计划制造两种产品I、II,下表为其耗费设备工时以及获利情况。问该公司制造两种家电各多少件能使利润最大。

项目 I II 每天可用能力
设备A/h 0 5 15
设备B/h 6 2 24
调试工序/h 1 1 5
利润/元 2 1

根据问题建立数学模型,并用单纯形法求得最终单纯形表:

cj 2 1 0 0 0
CB XB b x1 x2 x3 x4 x5
0 x3 15/2 0 0 1 5/4 -15/2
2 x1 7/2 1 0 0 1/4 -1/2
1 x2 3/2 0 1 0 -1/4 3/2
σj=cjZj 0 0 0 -1/4 -1/2

此时问题最优解为

X=(72,32,152,0,0)T

最优生产计划为生产72件产品I, 32件产品II。

例1:当产品I的利润下降至1.5元/件,产品II的利润增加至2元/件时,公司最优生产计划有何变化。

  • 将参数改变直接反映到最终单纯形表
    产品利润c1,c2,检验数
cj 1.5 2 0 0 0
CB XB b x1 x2 x3 x4 x5
0 x3 15/2 0 0 1 5/4 -15/2
1.5 x1 7/2 1 0 0 1/4 -1/2
2 x2 3/2 0 1 0 -1/4 3/2
σj=cjZj 0 0 0 1/8 -9/4
  • 检查原问题是否仍为可行解
    x1,x2,x3取值为正数,系数矩阵构成单位矩阵,基变量检验数全为0, 原问题是可行解。

  • 检查对偶问题是否仍为可行解
    x4的检验数大于零,对偶问题为非可行解。

  • 根据表中所列情况得出结论或决定继续计算的步骤
    由于原问题为可行解,对偶问题为非可行解,所以用单纯形法继续迭代求最优解。

cj 1.5 2 0 0 0
CB XB b x1 x2 x3 x4 x5
0 x4 6 0 0 4/5 1 -6
1.5 x1 2 1 0 -1/5 0
2 x2 3 0 1 1/5 0
σj=cjZj 0 0 -1/10 0 -3/2

此时生产计划调整为生产2件产品I、3件产品II。

1.2 约束条件右侧常数bi的灵敏度分析

约束条件的右侧常数表示资源的可利用程度,通过改变约束条件右侧常数可以分析资源的利用程度对决策变量解的影响。当约束条件右侧常数发生较大变动时,可能会改变最优解的取值范围。
当约束右侧常数项br有改变量Δbr时,不影响检验数的值,但是会影响最优解。将改变量后的右侧资源向量记为b,令 λ=CBB1=[λ1,λ2,,λm]表示右侧资源的影子价格。
XB=B1b0,则最优基不变,最优解的改变量ΔS=Δbrλr
若不满足XB=B1b0,则B不再是可行基,而变为正则基,因此将XB=B1bS=CBB1b代入原最优表,利用对偶单纯形法进行迭代,得到新的最优解。

例2:若原题中设备A和调试工序的每天能力不变,而设备B每天的能力增加到32h,分析公司最优计划的变化。

  • 将参数改变直接反映到最终单纯形表
    由题知,约束条件改变后,

b=[15245]Δb=[022]b=b+Δb

得到最终单纯形表中b列的变化,将其反映到最终单纯形表中。

cj 2 1 0 0 0
CB XB b x1 x2 x3 x4 x5
0 x3 35/2 0 0 1 5/4 -15/2
2 x1 11/2 1 0 0 1/4 -1/ 2
1 x2 -1/2 0 1 0 -1/ 4 3/2
σj=cjZj 0 0 0 -1/4 -1/2
  • 检查原问题是否仍为可行解
    x2取值变为-1/2, 原问题为非可行解
  • 检查对偶问题是否仍为可行解
    所有σj0,基变量检验数全为 0, 基变量的系数向量构成单位矩阵, 对偶问题为可行解
  • 按照表中所列情况得出结论或继续计算的步骤
    原问题为非可行解,对偶问题为可行解,因此使用对偶单纯形法进行迭代计算。
cj 2 1 0 0 0
CB XB b x1 x2 x3 x4 x5
0 x3 15 0 5 1 0 0
2 x1 5 1 1 0 0 1
0 x4 2 0 -4 0 1 -6
σj=cjZj 0 -1 0 0 -2

此时生产计划调整为生产5件产品I。

1.3 决策变量xj的灵敏度分析

决策变量的灵敏度分析可以评估决策变量值的改变对目标函数值和约束条件的违背程度产生的影响。通过改变决策变量的取值范围,可以判断最优解的稳定性和鲁棒性。
增加一个变量在实际问题中反映为增加一种新的产品。其分析步骤如下:

  • 计算σj=cjzj=cji=1maijyi
  • 计算Pj=B1Pj
  • σj0,原最优解不变,只需将计算得到的Pjσj直接写入最终单纯形表中;若σj>0,则按单纯形法继续迭代计算找出最优。

例3:设例题中的公司计划推出产品III,生产一件所需设备 A、设备 B及调试工序的时间分别为 3h、4h、2h,该产品的预期赢利为 3元/件,问该产品是否值得投产;如投产,对公司的最优生产计划有何变化。

设该公司x6件产品III,有c6=3,P6=(3,4,2)T

  • 将参数改变直接反映到最终单纯形表
    根据上述步骤计算新增变量在最终单纯形表中的检验数和系数列向量

σ3=3(0,14,12)[342]=1

P6=[15/415/201/41/201/43/2][342]=[702]

cj 2 1 0 0 0 3
CB XB b x1 x2 x3 x4 x5 x6
0 x3 15/2 0 0 1 5/4 -15/2 -7
2 x1 7/2 1 0 0 1/4 -1/ 2 0
1 x2 3/2 0 1 0 -1/4 3/2 2
σj=cjZj 0 0 0 -1/4 -1/2 1
  • 检查原问题是否仍为可行解
    基变量取值为正数,系数矩阵构成单位矩阵,基变量检验数全为 0, 原问题为可行解
  • 检查对偶问题是否仍为可行解
    σ6>0,对偶问题为非可行解
  • 按照表中所列情况得出结论或继续计算的步骤

原问题为可行解,对偶问题为非可行解,用单纯形法继续迭代计算

cj 2 1 0 0 0 3
CB XB b x1 x2 x3 x4 x5 x6
0 x3 51/4 0 7/2 1 3/8 -9/4 0
2 x1 7/2 1 0 0 1/4 -1/2 0
3 x6 3/4 0 1/2 0 -1/8 3/4 1
σj=cjZj 0 -1/2 0 -1/8 0 -5/4

由表可得,新产品值得投产,最优生产计划变为每天生产7/2件产品I, 51/4件产品 III。

二、灵敏度分析的步骤

灵敏度分析是研究与分析一个系统(或模型)的状态或输出变化对系统参数或周围条件变化的敏感程度的方法。‌它几乎在所有的运筹学方法以及在对各种方案进行评价时都是很重要的。灵敏度分析主要解决以下问题:

‌判断最优解的稳定性‌:当原始数据不准确或发生变化时,通过灵敏度分析来研究最优解的稳定性,决定哪些参数对系统或模型有较大的影响。
优化决策‌:通过分析各系数和约束条件的变化对模型的影响,决定是否需要调整决策。
方案评价‌:在有多种方案可供选择的情况下,通过灵敏度分析来评价不同方案的风险和可行性,选择最优方案。

灵敏度分析的步骤如下:

  • 将系数的改变通过计算反映到最终单纯形表
    按下列公式计算出由参数的变化而引起的最终单纯形表上有关数字的变化:

Δb=B1ΔbΔP=B1ΔPj

(cjzj)=cji=1maijyi

  • 检查原问题是否仍为可行解
    基变量的值中无负数、基变量的系数向量构成单位矩阵、基变量的检验数全为0,满足全部条件则为可行解。
  • 检查对偶问题是否仍为可行解
    所有σj0,基变量的检验数全为 0,基变量的系数向量构成单位矩阵,满足全部条件则对偶问题仍为可行解。
  • 按下表所列情况得出结论或决定继续计算的步骤
原问题 对偶问题 结论或继续计算的步骤
可行解 可行解 问题的最优解或最优基不变
可行解 非可行解 用单纯形法继续迭代求最优解
非可行解 可行解 用对偶单纯形法继续迭代求最优解
非可行解 非可行解 引进人工变量,编制新的单纯形表重新计算

三、灵敏度分析练习

例4: 试分析以下参数线性规划问题。当参数t0时的最优解变化。

maxz(t)=(3+2t)x1(5t)x2{x142x2123x1+2x218x1,x20

将此模型化为标准型

{maxz(t)=(3+2t)x1+(5t)x2x1+x3=42x2+x4=123x1+2x2+x5=18xj0,j=1,2,,5

t=0,用单纯形法求解,见下表3-1。

cj 3 5 0 0 0
CB XB b x1 x2 x3 x4 x5
0 x3 2 0 0 1 1/3 -1/3
5 x2 6 0 1 0 1/2 0
3 x1 2 1 0 0 -1/3 1/3
σj=cjZj 0 0 0 -3/2 -1

c的变化直接反映到最终表3-1中,得表3-2。

cj 3+2t 5t 0 0 0
CB XB b x1 x2 x3 x4 x5
0 x3 2 0 0 1 1/3 -1/3
5t x2 6 0 1 0 1/2 0
3+2t x1 2 1 0 0 -1/3 1/3
σj=cjZj 0 0 0 -3/2+7/6t -1+3/2t

t增大,t3/27/6=97时,首先出现σ0,在σ0,即 0t97时,得最优解(2,6,2,0,0)Tt=97为第一临界点。当t>97时,σ>0,这时x4作为换入变量。用单纯形法迭代一步,得表3-3。

cj 3+2t 5t 0 0 0
CB XB b x1 x2 x3 x4 x5
0 x4 6 0 0 3 1 -1
5t x2 3 0 1 -3/2 0 1/2
3+2t x1 4 1 0 1 0 0
σj=cjZj 0 0 9272t 0 52+12t

t继续增大t5/21/2=5时,首先出现σ0,在σ0,即97t5时,得最优解(4,3,0,6,0)Tt=5为第二临界点。当t>5时,σ>0,这时x3作为换入变量,用单纯形法迭代一步,得表3-4。

cj 3+2t 5t 0 0 0
CB XB b x1 x2 x3 x4 x5
0 x4 12 0 2 0 1 0
0 x5 6 0 2 -3 0 1
3+2t x1 4 1 0 1 0 0
σj=cjZj 0 5t 32t 0 0

表 3-4中t继续增大时,恒有σ<0,故当t5时,最优解为(4,0,0,12,6)T

四、灵敏度分析的Python实现

例5:由于分销和促销成本的差异,产品的利润也因销售渠道的不同而不同。此外,广告费用和人力成本也与销售渠道有关。下表给出了电子通信公司不同销售渠道的销售利润、广告费用和人工成本。广告预算5000美元,每个销售渠道的最大个人销售时间是1800小时,公司现阶段决定制造的产品数量600件,此外,全国连锁零售店要求最少销售150件产品。电子通信面临的问题是如何制定一个分销策略,使其总的销售利润最大。

分销渠道 单位已售产品的利润(美元) 单位已售产品的广告费用(美元) 单位已售产品的销售时间(小时)
航海器材经销店 90 10 2
商用器材经销店 84 8 2
全国连锁零售店 70 9 3
直接邮购 60 15

建立线性规划模型来解决电子通信公司如何制定分销策略使总的销售利润最大的问题。

  • 决策变量
    x1: 航海器材经销店的销售数量;x2: 商用器材经销店的销售数量;x3: 全国连锁零售店的销售数量;x4: 直接邮购的销售数量

  • 目标函数
    最大化总利润:Maximize Z=90x1+84x2+70x3+60x4

  • 约束条件

    • 广告费用约束:总广告预算为5000美元10x1+8x2+9x3+15x45000
    • 销售时间约束:每个分销渠道的最大销售时间为1800小时2x1+2x2+3x31800
    • 生产数量约束:总的生产数量为600件x1+x2+x3+x4=600
    • 最低销售数量约束:全国连锁零售店要求至少销售150件产品x3150
    • 非负约束:销售数量不能为负数x1,x2,x3,x40
  • 线性规划模型整合

MaxZ=90x1+84x2+70x3+60x4Subject to:10x1+8x2+9x3+15x450002x1+2x2+3x31800x1+x2+x3+x4=600x3150x1,x2,x3,x40

#上面模型求解和灵敏度分析
from ortools.linear_solver import pywraplp

# 创建求解器
solver = pywraplp.Solver.CreateSolver('GLOP')
if not solver:
    raise Exception("Solver not found.")

# 定义决策变量
x1 = solver.NumVar(0, solver.infinity(), 'x1')
x2 = solver.NumVar(0, solver.infinity(), 'x2')
x3 = solver.NumVar(0, solver.infinity(), 'x3')
x4 = solver.NumVar(0, solver.infinity(), 'x4')

# 目标函数
solver.Maximize(90 * x1 + 84 * x2 + 70 * x3 + 60 * x4)

# 约束条件
solver.Add(10 * x1 + 8 * x2 + 9 * x3 + 15 * x4 <= 5000)  # 广告费用约束
solver.Add(2 * x1 + 2 * x2 + 3 * x3 <= 1800)  # 销售时间约束
solver.Add(x1 + x2 + x3 + x4 == 600)  # 生产数量约束
solver.Add(x3 >= 150)  # 最低销售数量约束

# 求解
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL:
    print('Solution:')
    print(f'x1 = {x1.solution_value():.2f}')
    print(f'x2 = {x2.solution_value():.2f}')
    print(f'x3 = {x3.solution_value():.2f}')
    print(f'x4 = {x4.solution_value():.2f}')
    print(f'Total Profit = {solver.Objective().Value():.2f}')
    
    # 灵敏度分析
    print("\nSensitivity Analysis (Objective Coefficient Ranges):")
    x1_coef_lower = x1.ReducedCost() + 90
    x1_coef_upper = float('inf')  # 用无限表示无法具体量化的变化范围
    print(f"x1 coefficient range: [{x1_coef_lower:.2f}, {x1_coef_upper:.2f}]")
    
    x2_coef_lower = x2.ReducedCost() + 84
    x2_coef_upper = float('inf')
    print(f"x2 coefficient range: [{x2_coef_lower:.2f}, {x2_coef_upper:.2f}]")
    
    x3_coef_lower = x3.ReducedCost() + 70
    x3_coef_upper = float('inf')
    print(f"x3 coefficient range: [{x3_coef_lower:.2f}, {x3_coef_upper:.2f}]")
    
    x4_coef_lower = x4.ReducedCost() + 60
    x4_coef_upper = float('inf')
    print(f"x4 coefficient range: [{x4_coef_lower:.2f}, {x4_coef_upper:.2f}]")
else:
    print('The problem does not have an optimal solution.')
Solution:
x1 = 25.00
x2 = 425.00
x3 = 150.00
x4 = 0.00
Total Profit = 48450.00

Sensitivity Analysis (Objective Coefficient Ranges):
x1 coefficient range: [90.00, inf]
x2 coefficient range: [84.00, inf]
x3 coefficient range: [70.00, inf]
x4 coefficient range: [15.00, inf]

例6:已知线性规划模型及最优解

maxZ=x1+x2+3x3x1+x2+2x340x1+2x2+x320x2+x315x1,x2,x30

最优单纯形表

cj 1 1 3 0 0 0
CB XB b x1 x2 x3 x4 x5 x6
0 x4 5 0 -2 0 1 -1 -1
1 x1 5 1 1 0 0 1 -1
3 x3 15 0 1 1 0 0 1
σj 0 -3 0 0 -1 -2

b1的变化范围,使最优基不变;若b2增加10,求变化后的最优解。

from ortools.linear_solver import pywraplp

def main():
    # 创建求解器
    solver = pywraplp.Solver.CreateSolver('GLOP')  # 'GLOP' 是 Google 的线性求解器

    # 定义变量
    x1 = solver.NumVar(0, solver.infinity(), 'x1')
    x2 = solver.NumVar(0, solver.infinity(), 'x2')
    x3 = solver.NumVar(0, solver.infinity(), 'x3')

    # 定义目标函数
    solver.Maximize(x1 + x2 + 3 * x3)

    # 定义约束条件
    constraint_1 = solver.Add(x1 + x2 + 2 * x3 <= 40)
    constraint_2 = solver.Add(x1 + 2 * x2 + x3 <= 20)
    constraint_3 = solver.Add(x2 + x3 <= 15)

    # 求解问题
    status = solver.Solve()

    if status == pywraplp.Solver.OPTIMAL:
        print(f"最优解: x1 = {x1.solution_value()}, x2 = {x2.solution_value()}, x3 = {x3.solution_value()}")
        print(f"最优目标值: Z = {solver.Objective().Value()}")

        # 问题1:对约束1 (b1) 的灵敏度分析
        shadow_price_1 = constraint_1.dual_value()
        slack_1 = constraint_1.ub() - (x1.solution_value() + x2.solution_value() + 2 * x3.solution_value())
        
        # 允许变化范围:我们只能得到影子价格和松弛变量,范围需要手动推导
        if shadow_price_1 == 0:
            print(f"约束1的影子价格为 0,表明该约束不活跃,右侧 b1 的上限可以是无穷大")
            b1_lower_bound = 40 - slack_1
            print(f"约束1右侧常数 b1 的允许变化范围 (最优基不变): [{b1_lower_bound}, ∞)")
        else:
            b1_lower_bound = 40 - slack_1
            b1_upper_bound = 40 + slack_1 / abs(shadow_price_1)
            print(f"约束1的影子价格: {shadow_price_1}")
            print(f"约束1的松弛变量值: {slack_1}")
            print(f"约束1右侧常数 b1 的允许变化范围 (最优基不变): [{b1_lower_bound}, {b1_upper_bound}]")

        # 问题2:将约束2的右侧常数增加10,然后重新求解
        constraint_2.SetBounds(-solver.infinity(), 30)  # 右侧常数增加10
        solver.Solve()

        print(f"b2 增加10后的最优解: x1 = {x1.solution_value()}, x2 = {x2.solution_value()}, x3 = {x3.solution_value()}")
        print(f"b2 增加10后的最优目标值: Z = {solver.Objective().Value()}")
    else:
        print("未找到最优解")

if __name__ == '__main__':
    main()
最优解: x1 = 5.0, x2 = 0.0, x3 = 15.0
最优目标值: Z = 50.0
约束1的影子价格为 0,表明该约束不活跃,右侧 b1 的上限可以是无穷大
约束1右侧常数 b1 的允许变化范围 (最优基不变): [35.0, ∞)
b2 增加10后的最优解: x1 = 10.00, x2 = 0.0, x3 = 15.0
b2 增加10后的最优目标值: Z = 55.0

总结

线性规划的灵敏度分析是一种关键的工具,有助于决策者深入理解和有效应对参数变化所带来的影响,从而提高模型的实用性和可靠性。通过灵敏度分析,决策者可以评估目标函数系数、约束条件右端项以及约束系数矩阵元素的变化对最优解和目标值的影响。这种分析能够帮助在面对不确定性和环境变化时,保持决策的合理性和稳健性。在实际应用中,灵敏度分析的价值尤为显著。例如,在供应链管理中,当某些节点出现瓶颈或资源供给发生变化时,灵敏度分析可以帮助确定是否需要调整库存水平或改变供应链策略,从而避免资源浪费和生产中断。类似地,在投资组合优化中,灵敏度分析用于评估不同投资标的的收益系数变化对整体组合收益和风险的影响。这使得投资者能够在市场条件变化时及时调整投资策略,优化收益并控制风险。


参考文献

  1. 【线性规划(六)】对偶单纯形法与灵敏度分析
  2. 运筹说 第32期 | 对偶理论与灵敏度分析—灵敏度分析
  3. 对偶理论和灵敏度分析---参数线性规划
posted @   郝hai  阅读(658)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 全程不用写代码,我用AI程序员写了一个飞机大战
· DeepSeek 开源周回顾「GitHub 热点速览」
· 记一次.NET内存居高不下排查解决与启示
· 物流快递公司核心技术能力-地址解析分单基础技术分享
· .NET 10首个预览版发布:重大改进与新特性概览!
点击右上角即可分享
微信分享提示