pyomo建模学习笔记

1、建模要素

从运筹优化的角度，一般优化模型的构成包含三大要素：决策变量，约束，目标函数。在使用代码实现优化模型需要考虑：索引集合，参数集合，表达式。

2、例子：指派问题

$$\underset{x}{\min }\sum _{i\in I}\sum _{j\in J}{c}_{ij}{x}_{ij}$$
s.t.
$$\begin{gathered} \sum _{i\in I}{x}_{ij}=1,\forall j\in J \\ \sum _{j\in J}{x}_{ij}=1,\forall i\in I \\ {x}_{ij}\in [0,1],i\in I,j\in J \end{gathered}$$

• 索引集合：$I,J$
• 参数集合：$c_{ij},其中i\in I,j\in J$
• 决策变量： $x_{ij},其中i\in I,j\in J$
• 目标函数：${\sum }_{i\in I}{\sum }_{j\in J}{c}_{ij}{x}_{ij}$
• 约束条件：$$\begin{gathered} \sum _{i\in I}{x}_{ij}=1,\forall j\in J \\ \sum _{j\in J}{x}_{ij}=1,\forall i\in I \\ {x}_{ij}\in [0,1],i\in I,j\in J \end{gathered}$$

3、pyomo实现上面例子

3.1索引集合

3.1.1 Set(initialize,filter,validate)

• initialize:表示初始化参数，可以使用列表，元组和集合来给一个Set的初始化
• filter：表示一个过滤掉集合的条件，一般是一个函数表达式
• validate：表示对结合进行验证条件，一般是一个函数表达式

from pyomo.environ import *
model = ConcreteModel()
model.I = Set(initialize=[i for i in range(6)]) # 采用列表初始化集合
model.J = Set(initialize=('black','red','blue')) # 采用元组初始化集合
model.N = Set([i for i in range(3)])#等价于用initialize
model.M = Set([(i,j) for i in range(3) for j in range(3)])
def even_fun(model):
    return (i for model.I if i%2==0)
model.K = Set(initialize=even_fun)#可以通过一个函数来定义initialize

采用filter参数可以对集合的元素进行刷选，使用下面例子通过定义规则来实现model.L集合中只含有集合model.N没有的元素

def filter_rule(model,x):
    rerurn x not in model.N
model.L = Set(initialize = [i for i in range(6)],filter=filter_rule)

3.1.2 RangeSet(start,end,slice)

如果只需要纯数字作为索引的元素可以采用RangeSet会更加直接和方便，Set功能可以覆盖带掉RangeSet。使用方法如下：

model.C= RangeSet(1,10,2)#[1,3,5,7,9]
model.D = RangeSet(10)#[0,1,2,3,4,5,6,7,8,9,10]
model.E = RangeSet(1,10,4)#[1,5,9]
model.F = RangeSet(1,10)#[1,2,3,4,5,6,7,8,9,10]

3.1.3 集合的运算

model.I = [i for i in range(6)]
model.J = [i for i in range(6) if i%2 ==0]
model.K = [1,2]
model.E = model.I & model.J#交集（0,2,4）
model.F = model.I-model.J#差集（1,3,5）
model.G = model.I | model.J # 并集（0，1,2,3,4,5）
model.H = model.J * model.K #笛卡尔积（((0, 1), (0, 2), (2, 1), (2, 2), (4, 1), (4, 2))）

3.2 参数集合

Param(initialize,default,validate)

• initialize:表示初始化参数，可以采用列表，元组或字典来做参数初始化
• default：表示初始化参数的默认值，initialize没有覆盖的索引default才会起作用
• validate：表示对集合进行验证条件，一般是一个函数表达式

3.2.1 简单的赋值法初始化参数

from pyomo.envirn import *
model = ConcreteModel()
model.I = Set(initialize = [i for i in range(4)])
model.J = Set(initialize = [i for i in range(4)])
model.c = Param(model.I,model.J,default = 1)#参数c是一个4*4的矩阵并且初始值全部为1

3.2.2 字典初始化参数

c_data = {(0,0):1,(1,2):33,(3,3):20}# 采用字典初始化 表示 c[0,0]=1,c[1,2]=33,c[3,3]=20
model.c1 = Param(model.I,model.J,initaialize=c_data,default=0)#default = 0 表示其余元素都为 0

3.2.3 函数初始化参数

采用函数初始化参数相比前两种方法可以非常的灵活

def c_init(model,i,j):
    if i==j:
        return i*i
    else:
        return 0
model.c2 = Param(model.I,model.J,initialize = c_init)

3.2.4 validate 验证参数取值正确性

借用 c_validate 函数可以方便验证参数选取是否满足一些硬性条件，如果不满足条件则会直接报错提示。这是一个非常好的预防bug提升代码鲁棒性的功能，同时validate并不影响我们对数学模型的实现。

c_data = {(1,1): 4, (1,2): 33, (3,3): 20}
def c_validate(model, v, i, j): # 保证所有c_{ij}的值都要大于3
    return v > 3.

model.c3 = Param(model.I, model.J, initialize=c_data, validate=c_validate)

在上面函数中，如果$c_{ij}$的值小于3，validate函数会报错提示我们

3.2.5 参数常用操作

print(len(model.c))#9 参数元素的个数
print(1 in model.c)#false
print((0,0) in model.c)#True
print([v for v in model.c])# [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)] 遍历集合中的元素

3.3 决策变量

Var(index_set,within,bounds,initialize,domain)

• index_set:表示决策变量的索引
• within：表示决策变量的类型
• bounds：表示决策变量上下界
• initialize：表示决策变量的初值
• domain 用于指定变量的数学类型，实数(Reals)、整数(Integers)和布尔值(Boolean)等

3.3.1 一维决策变量的定义

from pyomo.environ import *
model = ConcreteModel()
model.x = Var(within=Reals,bounds=(0,6),initialize=2)#表示实数变量，上界是6，下界为0，初值设为2

3.3.2 多维决策变量定义

def bounds_fun(model,i,j):
    return (0,1)#上界为1，下界为0
model.I = Set(initialize=[i for i in range(4)])
model.J = Set(initialize = [j for j in range(4)])
model.y = Var(model.I,model.J,within=PositiveIntegers,bounds = bounds_fun)#表示决策变量y维数是4*4 正整数变量，上下界为1,0
def bouns_fun1(model,i,j):
    return (10,None)#没有上界
model.z = Var(model.I,model.J,within=PositiveIntegers,bounds = bounds_fun1)

within 可以选择的类型：

• Reals： 实数
• PositiveReals：正实数
• NonPositiveReals：非正实数
• NegativeReals：负实数
• NonNegativeReals：非负实数
• PercentFraction：0-1之间实数 = UnitInterval
• Integers：整数
• PositiveIntegers：正整数
• NonPositiveIntegers：非正整数
• NegativeIntegers:负整数
• NonNegativeIntegers:非负整数
• Boolean: 布尔变量
• Binary: 0,1变量

3.3.3 决策变量赋初值

在pyomo中通过initialize参数来实现给决策变量赋初值的功能

model.x = Var(model.I,model.J,initialize={(0,0):1.2,(1,2):30,(2,2):21}) # 字典赋初值给决策变量
model.y = Var(model.I,model.J,initialize=2.0)#初值全是2.0
def g_init(model, i, j):
    return i+j

model.x_init3 = Var(model.I, model.J, initialize = g_init) # 通过函数g_init 给定初值
print(value(model.x_init1[0,0])) # 1.2
print(value(model.x_init1[1,2])) # 30
print(value(model.x_init2[0,0])) # 2.0
print(value(model.x_init2[3,3])) # 2.0
print(value(model.x_init3[3,3])) # 6.0

3.3.4 决策变量设置固定值

为了降低问题的规模和难度，经常采用固定一部分变量，只优化部分变量，在pyomo中提供了固定变量的方法，一旦变量被固定之后在求解过程中该变量视为常数。

model.x.fix(1)#表示变量x固定为1
print(model.x[0,0].fixed)#True 表示该变量是固定变量
print(value(model.x))#1

3.4 目标函数

Objective(index_set,obj_rule,sense)

• index_set:表示目标函数的索引集合
• obj_rule：表示目标函数定义规则
• sense：maximize表示最大化问题，minimize表示最小化问题

3.4.1 定义目标函数

from pyomo.environ import *
model = ConcreteModel()
model.I = Set(initialize=[i for i in range(4)])
model.J = Set(initialize=[j for j in range(4)])
model.x = Var(model.I,model.J,domain=Reals)
model.obj = Objective(rule = 2*model.x[1,1]+3*model.x[1,2],sense=minimize)

summation 函数是pyomo中可以实现两个向量快速求内积的操作，
$$\sum _{i\in I}\sum _{j\in J}{p}_{ij}{x}_{ij}$$

model.p=Param(model.I,model.J,default=2)
def obj_rule(model):
    return summation(model.p,model.x)
model.obj1 = Objective(rule=obj_rule,sense=maximize)

Pyomo 中可以建模各种各样的目标函数，如下的目标函数中包括了三角函数，指数函数和对数函数。

def obj_rule(model):
    expr = 0
    for i in model.I:
        for j in model.J:
            expr += sin(model.x[i,j]) + model.x[i,j]**2 + 2**(model.x[i,j])+log(model.x[i,j])
    return expr
model.obj2= Objective(rule=obj_rule,sense=maximize)

带有索引的目标函数的定义

def e_rule(model,i,j):
    return model.x[i,j]**2
model.e1=Objective(model.I,model.J,rule=e_rule)#对于任意一个i，j求x[i,j]的平方

使用Skip实现目标函数的灵活定义，如果某些索引要跳过，可以使用Objective.Skip,

def e_rule(model,i,j):
    if i==j:
        return Objective.Skip
    else:
        return model.x[i,j]**2
model.e2 = Objective(model.I,model.J,rule=e_rule,sense=minimize)

打印目标函数

print(model.obj1.expr)#2*x[1,1]+3*x[2,2] 返回目标函数的表达式
print(model.obj1.sense)#-1表示极大化，1表示极小化
print(value(model.obj1))#返回目标函数的值

3.5 约束

Constraint(index_set,expr,rule)

• index_set:表示约束索引集合
• rule : 表示约束定义规则
• expr ： 约束表达式

from pyomo.environ import *
model = ConcreteModel()
model.I = Set(initialize = [i for i in range(3)])
model.J = Set(initialize = [j for j in range(3)])
model.x = Var(model.I,model.J,initialize=1,within=Reals)
model.c1 = Constraint(expr = model.x[0,0] + 2*model.x[0,3]>=2)#不等式约束
model.c2 = Constraint(expr=3*model.x[0,0]+5*model.x[1,1]==10)#等式约束
print(value(model.c1.body))#3返回c1的值
print(model.c1.lslack())#1 返回c1约束与下限的值（lower bound slack variable）
print(model.c1.uslack())#+inf 返回c1约束与上限的差（upper bound slack variable）
print(model.c2.lslack()) # -2.0 
print(model.c2.uslack()) # 2.0 

借助函数定义约束，如${\sum }_{j\in J}{p}_{ij}{x}_{ij}\le 2,\forall i\in I$

def c_rule(model,i):
    if i==0:
        return Constraint.Skip
    else:
        return sum([model.p[i,j] * model.x[i,j] for j in model.J])<=2
model.c3 = Constraint(model.I,rule=c_rule)

激活约束和停用约束

model.c2.deactivate()#c2约束失效
model.c2.activate()#c2约束激活

调用 deactivate() 后可以让约束c2不起作用了，如果想再添加回c2约束，再调用activate()即可

3.6 表达式

Expression(expr,rule)

• expr:表达式
• rule:定义表达式的函数

from pyomo.environ import *
model = ConcreteModel()
model.I = Set(initialize = [i for i in range(3)])
model.x = Var(model.I, within = Reals, initialize = 1.0)
model.e = Expression(expr = model.x[0] + 3*model.x[1]) # 建立表达式e=x[0]+3x[1]
print(value(model.e)) # 4.0
model.e.set_value(model.x[1]) # 重新给表达式赋值 e = x[1]
print(value(model.e)) # 1.0
#定义多个表达式，xi*xi + exp(xi)
def gen_expr_rule(model,i):
    return model.x[i]* model.x[i] + exp(model.x[i])

model.e2 = Expression(model.I, rule=gen_expr_rule)

4 求解

Pyomo 只能完成建模的任务，求解的还是需要将 Pyomo 的模型导入到专业的优化求解器来实现。如下代码就是调用 Gurobi 求解其来实现求解：

from pyomo.environ import *
model = ConcreteModel()
model.I = Set([i for i in range(4)])
model.J= Set([i for i in range(4)])
model.x = Var(model.I, model.J, within = Reals, bounds = (0,1))
model.c = Param(model.I, model.J, initialize = {(1,1):2, (2,2):3}, default = 1)
def obj_rule(model):
    return summation(model.c,model.x)
model.obj = Objective(rule=obj_rule,sense=minimize)
def con_rule1(model,i):
    return sum([model.x[i,j] for j in model.J]) ==1
model.con1 = Constraint(model.I,rule= con_rule1)
def con_rule2(model,j):
    return sum([model.x[i,j] for i in model.I]) ==1
model.con2 = Constraint(model.J,rule= con_rule2)
model.dual = Suffix(direction=Suffix.IMPORT)#对偶变量定义
opt = SolverFactory('gurobi')#指定gurobi作为求解器
# opt = SolverFactory('cplex')#指定cplex作为求解器
# opt = SolverFactory('scip')#指定scip作为求解器
solution = opt.solve(model)#调用求解器求解

注意在调用求解器之前要确保求解器已经成功安装。Pyomo 目前支持的求解器种类也有很多种主要包括：Gurobi，Cplex，SCIP, CBC, Ipopt, BARON, Moske, Xpress, GLPK等等。调用方式和上面例子中的Gurobi 调用方式类似，只需要替换相关求解器的名字即可。

打印结果

solution.write()#写入求解结果
x_opt = np.array([value(model.x[i,j]) for i in model.I for j in model.J]).reshape((len(model.I),len(model.J)))
obj_vale = value(model.obj)#提取目标函数的值
print("optimum point: \n {} ".format(x_opt))
print("optimal objective: {}".format(obj_values))

如果需要使用对偶变量，需要在求解前就要定义好对偶变量，如
model.dual = Suffix(direction=Suffix.IMPORT)这样在求解完成之后，可以使用羡慕方法获得对偶变量的值

for c in model.con1:
    print(model.dual[model.con1[c]])#遍历约束con1 的对偶变量

print(solution.solver.termination_condition) # 终止条件 一般包括三种 optimal, feasible, infeasible
print(solution.solver.termination_message) # 求解得到的信息 一般为一行字符串信息
print(solution.solver.status) # 表示求解的状态 包括 ok, warning, error, aborted, or unknown

5 例子

$$\underset{x}{\min }\sum _{i\in I}\sum _{j\in J}{c}_{ij}{x}_{ij}$$
s.t.
$$\begin{gathered} \sum _{i\in I}{x}_{ij}=1,\forall j\in J \\ \sum _{j\in J}{x}_{ij}=1,\forall i\in I \\ {x}_{ij}\in [0,1],i\in I,j\in J \end{gathered}$$
完整代码如下：

from pyomo.environ import *
model = ConcreteModel()
model.I = Set(initialize=[i for i in range(4)])
model.J = Set(initialize=[j for j in range(4)])
model.x = Var(model.I,model.J,within=Reals,bounds=(0,1))
model.c = Param(model.I, model.J, initialize = {(1,1):2, (2,2):3}, default = 0)
def obj_rule(model):
    return summation(model.c,model.x)
model.obj = Objective(rule=obj_rule,sense=minimize)
def con_rule1(model,i):
    return sum([model.x[i,j] for j in model.J]) ==1
model.con1 = Constraint(model.I,rule= con_rule1)
def con_rule2(model,j):
    return sum([model.x[i,j] for i in model.I]) ==1
model.con2 = Constraint(model.J,rule= con_rule2)
opt =SolverFactory('gurobi')
solution = opt.solve(model)
solution.write()
x_opt = np.array([value(model.x[i,j]) for i in model.I for j in model.J]).reshape((len(model.I), len(model.J))) # 提取最优解
obj_values = value(model.obj) # 提取最优目标函数值
print("optimum point: \n {} ".format(x_opt))
print("optimal objective: {}".format(obj_values))

参考

https://pyomo.readthedocs.io/en/stable/