python 线性规划_鲁棒线性规划问题的建模及求解

最新推荐文章于 2024-08-14 21:31:57 发布

weixin_39723920

最新推荐文章于 2024-08-14 21:31:57 发布

阅读量593

点赞数

文章标签： python 线性规划

参考资料

留德华叫兽：优化 | 鲁棒优化基础zhuanlan.zhihu.com

https://www2.isye.gatech.edu/~nemirovs/MP_RO_2002.pdfwww2.isye.gatech.edu

鲁棒优化(Robust optimization)

鲁棒优化是指，定义问题的数据不正确或者不确定的情况下，可以返回信赖的结果的优化问题的建模技术及算法。

由于优化的结果对于外部的扰乱表现的非常不稳定，存在完全无法信赖的结果的情况。比如(Aharon Ben-Tal and Arkadi Nemirovski )，非常有名的线性规划问题的 benchmark (NETLIB)中，对定义的线性约束的系数做比较小的扰动的情况下，问题可能变得不可行。

二次锥优化问题(Second-order cone optimization problem)

这类优化问题，如他的名字一样，因为约束里有二次锥约束。

最基本的3次元的二次锥为

3次元的二次锥

看它长得像什么呢？冰激凌，所以叫它冰激凌锥。别急它还有个更加高大上的名字，叫做洛伦兹锥(Lorentz cone)。锥的种类很多，大家可详细搜索一下。

混合问题(Product mix problem)的建模

这里我们考虑一个古典的线性规划问题-混合问题，它可描述为:

这里有一个工厂。它要把四种原料混合之后，制造一种产品。对于原料包含三种成分。我们想混合成为，使成分1百分之20以上，成分2百分之30以上，成分3百分之20以上。各原料成分的含有比例如图，原料单价每吨分别是5,6,8,20元。那么，怎么混合原料使得用最小费用来制造产品？

各原料成分含有比例

我们可把混合问题用线性规划问题进行建模如下:

中，原料

的价格是

,在原料

里包含成分

的比率为

,产品中应该包含的成分

的比率的下限为

。表示原料

的混合比率的变量为

。这样用求解器直接求解即可。

但是，实际中，原料里包含的成分的比率能够正确的确定的情况很少，多少会伴随有点误差。

在这里我们考虑原料

里包含成分

的比率

伴随有误差

。另外，我们需要对这个误差在一定范围内定义。这里我们考虑

表示限制误差量的参数。

那么，考虑上述的误差，

中的约束条件

可重新整理为

进而可得

或

这里，在式子左边，取

的地方是表示在4次元的半径为

的球上的线性函数

的最小化。但是，变量是

。

什么？什么？4次元!

为了问题的简单化，我们考虑2次元的半径为1的球(圆)上的线性函数

最小化的情况。

从上图可知，在圆上使

最小的点

有

进而推广到现在这个问题上，可得线性函数在球上的最优解为

但是，

。所以，我们的约束进而可表示为

最终，

中含有误差的约束可表示为二次锥约束，如下

Reformulation后，鲁棒优化问题可用二次锥最优化问题建模为

如果这里

那么

。

的求解

为二次锥优化问题，Gurobi (内点法已经OK了，From Version 5.0)是可以 handle。这里，我们用编程语言Python调用Gurobi的API进行求解。

#code reference: Prof. Kubo
from gurobipy import *

def prodmix(I,K,a,p,epsilon,LB):
    """prodmix:  robust production planning using soco
    Parameters:
        I - set of materials
        K - set of components
        a[i][k] -  coef. matrix
        p[i] - price of material i
        LB[k] - amount needed for k
    Returns a model, ready to be solved.
    """

    model = Model("robust product mix")
    x,rhs = {},{}
    for i in I:
        x[i] = model.addVar(vtype="C", name="x(%s)"%i)
    for k in K:
        rhs[k] = model.addVar(vtype="C", name="rhs(%s)"%k)
    model.update()

    model.addConstr(quicksum(x[i] for i in I) == 1)
    for k in K:
        model.addConstr(rhs[k] == -LB[k]+ quicksum(a[i,k]*x[i] for i in I) )
        model.addConstr(quicksum(epsilon*epsilon*x[i]*x[i] for i in I) <= rhs[k]*rhs[k])

    model.setObjective(quicksum(p[i]*x[i] for i in I), GRB.MINIMIZE)

    model.update()
    model.__data = x,rhs
    return model


def make_data():
    a = { (1,1):.25, (1,2):.15, (1,3):.2,
          (2,1):.3,  (2,2):.3,  (2,3):.1,
          (3,1):.15, (3,2):.65, (3,3):.05,
          (4,1):.1,  (4,2):.05,  (4,3):.8
          }
    epsilon = 0.01#0.02
    I,p = multidict({1:5, 2:6, 3:8, 4:20})
    K,LB = multidict({1:.2, 2:.3, 3:.2})
    return I,K,a,p,epsilon,LB


I,K,a,p,epsilon,LB  = make_data()
model = prodmix(I,K,a,p,epsilon,LB)
model.optimize()
print ("obj:",model.ObjVal)
x,rhs = model.__data
for i in I:
    print (i,x[i].X)

如果感兴趣的可以做一下当

和

的最优解以及目标函数值的差别。

如果有什么疑问或者建议，可以给我发邮件。➡ zhaoyou728 atoutlook.com