我们之前都是每个决策变量比较少，可以很方便的写出来（比如只有x1,x2两个决策变量），但如果变量比较多（x1-xi，i是比较大的整数的情况下），也就是我们的目标函数和约束条件出现了连加符号∑的情况，我们的程序该怎么写呢？本篇以啤酒混合策略、寒冷天气混合燃料策略、饮食营养策略三个例子进行讲解。
一、实践1:啤酒混合策略(Jenchura, 2017)
1、问题描述
一家酿酒厂收到了100加仑4% ABV(酒精含量)啤酒的订单。啤酒厂现有啤酒A为4.5% ABV，每加仑成本为0.32美元，啤酒B为3.7% ABV，每加仑成本为0.25美元。水也可用作混合剂，每加仑成本为0.05美元。找到满足客户要求的最低成本组合。（可以有更多组成部分component）
首先是用一种通用的方式表示问题数据，如果需要，可以扩展这种方式来包含额外的混合组件。这里我们使用一个材料字典，每个键表示一个混合剂。对于每个键，都有一个包含每个混合组件属性的子字典。
data = {
'A': {'abv': 0.045, 'cost': 0.32},
'B': {'abv': 0.037, 'cost': 0.25},
'W': {'abv': 0.000, 'cost': 0.05},
}2、数学建模
目标函数：
如果我们让下标c表示组成成分c（即component）集合中的一个混合分量，并将混合中使用的c的体积表示为xc，则混合的代价为：
其中Pc为c的单位体积价格。使用上面定义的Python数据字典，Pc的价格由data[c]['cost']给出。比如data['A']['cost']就是A的成本。
体积约束：
顾客的要求是生产出总的体积v。假设有理想解，约束条件为：
Xc是组成部分的体积。
产品组成成分的约束
产品的组成按体积规定为酒精的4%，记为A上横线，Ac是组成成分c的酒精浓度。这个约束可以这样表示：
把这个约束改写成线性约束：
3、编程实现
安装cbc求解器，也可以使用之前使用的glbk求解器:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import shutil
import sys
import os.path
if not shutil.which("pyomo"):
!pip install -q pyomo
assert(shutil.which("pyomo"))
if not (shutil.which("cbc") or os.path.isfile("cbc")):
if "google.colab" in sys.modules:
!apt-get install -y -qq coinor-cbc
else:
try:
!conda install -c conda-forge coincbc
except:
pass
assert(shutil.which("cbc") or os.path.isfile("cbc"))模型实现：
import pyomo.environ as pyomo
vol = 100
abv = 0.040
def beer_blend(vol, abv, data):
C = data.keys()
model = pyomo.ConcreteModel()
model.x = pyomo.Var(C, domain=pyomo.NonNegativeReals)
model.cost = pyomo.Objective(expr = sum(model.x[c]*data[c]['cost'] for c in C)) # 目标函数是成本的连加
model.vol = pyomo.Constraint(expr = vol == sum(model.x[c] for c in C))
model.abv = pyomo.Constraint(expr = 0 == sum(model.x[c]*(data[c]['abv'] - abv) for c in C))
solver = pyomo.SolverFactory('cbc') # 也可以采用glbk求解器
solver.solve(model)
print('Optimal Blend')
for c in data.keys():
print(' ', c, ':', model.x[c](), 'gallons')
print()
print('Volume = ', model.vol(), 'gallons')
print('Cost = $', model.cost())
beer_blend(vol, abv, data)
二、实践2:寒冷天气燃料设计问题
1、问题描述
酒精炉的燃料通常是变性酒精denatured alcohol。变性酒精，有时也被称为甲基化酒精，通常是乙醇和其他酒精及化合物的混合物，旨在使其不适合人类消费。一个制造商产品的MSDS描述描述了乙醇和甲醇的大约50%的混合。
酒精炉的问题是在寒冷的天气下很难点燃。这个问题的目的是设计一种可替代的寒冷天气燃料（设计一种新燃料可以替代变性酒精，低温下（-10度）还能燃烧，需要此时蒸汽压足够大；高温下（30度）不会爆炸，此时需要蒸汽压不超过变性酒精在30度的蒸汽压），它可以混合其他材料，通常可以从硬件或家装店获得。
计算蒸汽压公式
以下数据收集的潜在燃料。数据包括价格($/加仑)和参数，用于预测Antoine蒸汽压方程，下标s表示物质，温度T的单位是摄氏度，压力P的单位是mmHg。需要额外的分子量和比重信息来表示体积分数的最终结果。物质s主要参与计算的性质是A/B/C，MW和SG主要是根据摩尔分数（也就是这个问题的决策变量xs）计算质量分数、体积分数时会用到。
data = {
'ethanol' : {'MW': 46.07, 'SG': 0.791, 'A': 8.04494, 'B': 1554.3, 'C': 222.65},
'methanol' : {'MW': 32.04, 'SG': 0.791, 'A': 7.89750, 'B': 1474.08, 'C': 229.13},
'isopropyl alcohol': {'MW': 60.10, 'SG': 0.785, 'A': 8.11778, 'B': 1580.92, 'C': 219.61},
'acetone' : {'MW': 58.08, 'SG': 0.787, 'A': 7.02447, 'B': 1161.0, 'C': 224.0},
'xylene' : {'MW': 106.16, 'SG': 0.870, 'A': 6.99052, 'B': 1453.43, 'C': 215.31},
'toluene' : {'MW': 92.14, 'SG': 0.865, 'A': 6.95464, 'B': 1344.8, 'C': 219.48},
}第一步是在一个典型的操作温度范围内测定变性酒精的蒸汽压。为此，我们假设变性酒精是乙醇和甲醇的40/60(摩尔分数)的混合物。
def Pvap(T, s):
return 10**(data[s]['A'] - data[s]['B']/(T + data[s]['C'])) # 即上述的Antoine蒸汽压方程
def Pvap_denatured(T):
return 0.4*Pvap(T, 'ethanol') + 0.6*Pvap(T, 'methanol') # 假设变性酒精是乙醇和甲醇的40/60(摩尔分数)的混合物，那么它的蒸汽压是多少，画出下图（它蒸汽压和温度T的关系）
T = np.linspace(0, 40, 200) # 温度范围0-40摄氏度，均匀划分成200个点
plt.plot(T, Pvap_denatured(T))
plt.title('Vapor Pressure of denatured alcohol')
plt.xlabel('temperature / °C')
plt.ylabel('pressure / mmHg')
print("Vapor Pressure at 0C =", round(Pvap_denatured(0),1), "mmHg")
plt.grid(True)
产品要求
我们寻求一种能在0°C或更低的温度下增加蒸汽压力的低温燃料，并能在更高的操作温度下保证酒精炉的安全正常运行。
为此，我们寻求一种蒸汽压至少为22毫米汞柱、温度最低、不大于30°C及以上的变性酒精蒸汽压的常用液体混合物。
for s in data.keys():
plt.plot(T, Pvap(T,s))
plt.plot(T, Pvap_denatured(T), 'k', lw=3)
plt.legend(list(data.keys()) + ['denatured alcohol'])
plt.title('Vapor Pressure of selected compounds')
plt.xlabel('temperature / °C')
plt.ylabel('pressure / mmHg')
plt.grid(True)
2、数学模型
优化模型是创建一个混合物，在-10°C时蒸汽压最大，而在30°C及以上时蒸汽压小于或等于变性酒精。
优化模型中的决策变量与xs对应，xs是可用物种s集合中每个物种s∈s的摩尔分数（mole fractions）。根据定义，摩尔分数必须满足
其目的是在低温条件下，即-10°C时，使蒸汽压达到最大（使之在低温环境下酒精灯还可以燃烧），而在30°C时，使蒸汽压小于或等于变性酒精。利用理想混合物的拉乌尔定律，
3、编程实现
m = pyomo.ConcreteModel()
S = data.keys()
m.x = pyomo.Var(S, domain=pyomo.NonNegativeReals)
def Pmix(T):
return sum(m.x[s]*Pvap(T,s) for s in S)
m.obj = pyomo.Objective(expr = Pmix(-10), sense=pyomo.maximize) # 目标函数
m.cons = pyomo.ConstraintList()
m.cons.add(sum(m.x[s] for s in S)==1) # 摩尔分数之和等于1
m.cons.add(Pmix(30) <= Pvap_denatured(30)) # 30°C时，混合燃料蒸汽压小于或等于变性酒精蒸汽压
m.cons.add(Pmix(40) <= Pvap_denatured(40)) # 40°C时
solver = pyomo.SolverFactory('cbc')
solver.solve(m)
print("Vapor Pressure at -10°C =", m.obj(), "mmHg")
T = np.linspace(-10,40,200)
plt.plot(T, Pvap_denatured(T), 'k', lw=3)
plt.plot(T, [Pmix(T)() for T in T], 'r', lw=3)
plt.legend(['denatured alcohol'] + ['cold weather blend'])
plt.title('Vapor Pressure of selected compounds')
plt.xlabel('temperature / °C')
plt.ylabel('pressure / mmHg')
plt.grid(True)也就是说，在满足3个约束条件下，混合燃料最大蒸汽压是17.48 mmHg , 混合燃料和变性酒精的对比如下:
s = data.keys()
results = pd.DataFrame.from_dict(data).T
for s in S:
results.loc[s,'mole fraction'] = m.x[s]()
MW = sum(m.x[s]()*data[s]['MW'] for s in S)
for s in S:
results.loc[s,'mass fraction'] = m.x[s]()*data[s]['MW']/MW
vol = sum(m.x[s]()*data[s]['MW']/data[s]['SG'] for s in S)
for s in S:
results.loc[s,'vol fraction'] = m.x[s]()*data[s]['MW']/data[s]['SG']/vol
results
混合燃料的各组成成分的摩尔分数、质量分数、体积分数如下：
三、实践3:饮食营养问题(Danzig，1990)
1、问题描述
饮食问题的目标是选择以最低成本满足日常营养需求的食物。这个问题可以表述为一个线性程序，它限制了饮食中卡路里的数量和维生素、矿物质、脂肪、钠和胆固醇的数量。Danzig(1990)指出，饮食问题的动机是美军希望在提供健康饮食的同时将战场上的士兵的饮食成本最小化。也就是每个士兵的最小饮食成本是多少？
2、数学模型
数学模型如下：
- F是食物集合，N是营养集合，aij表示i食物含有多少j营养，ci是每单位i食物的成本，xi是消耗了多少单位食物，这个问题就是要在满足营养需求的情况下最小化成本指出（看目标函数objective）。
- Nmin/max j 是营养j的最小/最大需求，这个范围就是营养范围；
- Vi是单位食物i的体积，Vmax是一个士兵能携带的最大食物体积；
- 约束条件就是不大于最大体积，营养在营养范围之内。
3、编程实现
定义参数
from pyomo.environ import *
infinity = float('inf')
model = AbstractModel() # 这是一个抽象的Pyomo模型，因此需要提供set和参数值来初始化模型，见下面的“问题数据：通过dat文件读取”。
# Foods
model.F = Set()
# Nutrients
model.N = Set()
# Cost of each food
model.c = Param(model.F, within=PositiveReals) # ci
# Amount of nutrient in each food
model.a = Param(model.F, model.N, within=NonNegativeReals) # aij
# Lower and upper bound on each nutrient
model.Nmin = Param(model.N, within=NonNegativeReals, default=0.0) # Nmin-j
model.Nmax = Param(model.N, within=NonNegativeReals, default=infinity) # Nmax-j
# Volume per serving of food
model.V = Param(model.F, within=PositiveReals) # Vi
# Maximum volume of food consumed
model.Vmax = Param(within=PositiveReals) # Vmax这些参数声明中使用了within选项来定义参数的预期属性。此信息用于对用于初始化参数组件的数据执行错误检查。
决策变量
Var组件用于定义决策变量:
# Number of servings consumed of each food
model.x = Var(model.F, within=NonNegativeIntegers) # 非负整数
目标函数
# Minimize the cost of food that is consumed
def cost_rule(model):
return sum(model.c[i]*model.x[i] for i in model.F)
model.cost = Objective(rule=cost_rule)
约束条件
# Limit nutrient consumption for each nutrient
def nutrient_rule(model, j):
value = sum(model.a[i,j]*model.x[i] for i in model.F)
return model.Nmin[j] <= value <= model.Nmax[j]
model.nutrient_limit = Constraint(model.N, rule=nutrient_rule)
# Limit the volume of food consumed
def volume_rule(model):
return sum(model.V[i]*model.x[i] for i in model.F) <= model.Vmax
model.volume = Constraint(rule=volume_rule)以上代码记在diet.py文件里。
问题数据：通过dat文件读取
例子一我们使用两个字典表示问题数据，这里我们使用diet.dat文件读取：
求解模型
使用命令行求解：
!pyomo solve --solver=glpk diet.py diet.dat
获取结果文件：
!cat results.yml
这个解决方案显示，每天花费大约15美元，一个人可以吃4个芝士汉堡、5根薯条、1个鱼肉三明治和4罐牛奶。
参考
【3】pyomo例子3:https://nbviewer.jupyter.org/github/Pyomo/PyomoGallery/blob/master/diet/DietProblem.ipynb