Docplex入门（1）——线性规划
| 阿里云国内75折 回扣 微信号：monov8 |
| 阿里云国际，腾讯云国际，低至75折。AWS 93折 免费开户实名账号 代冲值 优惠多多 微信号：monov8 飞机：@monov6 |
Docplex入门——线性规划
定义及前置知识
- 线性规划解决最大/小化一个线性目标函数所有决策变量都是连续且要求目标和约束都由线性表达式组成
- 线性表达式linear_expression
∑ a i x i 亦 可 写 为 矩 阵 形 式 A X \sum{a_ix_i} 亦可写为矩阵形式AX ∑aixi亦可写为矩阵形式AX
其中要求A为常数向量已知X为自变量待求解且连续对比一下非线性的例子
1.两个/多个变量的乘法
2.指数
3.对数
4.绝对值
..
-
线性约束
l i n e a r − e x p r e s s i o n s y m b o l l i n e a r − e x p r e s s i o n linear_ -expression \quad symbol \quad linear_-expression linear−expressionsymbollinear−expression
其中symbol仅能取 = = =、 ≤ \le ≤、 ≥ \ge ≥ -
一般形式
min C x t s . t . A x ≥ B x ≥ 0 \min C^t _x \\s.t. \quad Ax \ge B \\x \ge 0 minCxts.t.Ax≥Bx≥0
其中x为大小为n的向量A为m*n的系数矩阵B是大小为m的常向量。
线性规划样例
-
问题描述
要生产台式和移动两类电话要求每种类型电话至少生产100台且最大化利润。但是公司生产能力有限制造过程分为装配和涂装2个步骤需要在不超过生产能力情况下计算每款手机最优生产数量。 -
建模
决策变量 A 、 B 各 自 的 数 量 A、B各自的数量 A、B各自的数量
目标 m a x i m i z e : 12 d e s k − p r o d u c t i o n + 20 c e l l − p r o d c t i o n maximize: 12desk_-production + 20cell_-prodction maximize:12desk−production+20cell−prodction
约束条件
d e s k − p r o d u c t i o n ≥ 100 c e l l − p r o d u c t i o n ≥ 100 0.2 d e s k − p r o d u c t i o n + 0.4 c e l l − p r o d u c t i o n ≤ 400 0.5 d e s k − p r o d u c t i o n + 0.4 c e l l − p r o d u c t i o n ≤ 490 \ desk_-production \ge 100 \\ cell_-production \ge 100 \\ 0.2desk_-production + 0.4cell_-production \le 400 \\ 0.5desk_-production + 0.4cell_-production \le 490 desk−production≥100cell−production≥1000.2desk−production+0.4cell−production≤4000.5desk−production+0.4cell−production≤490 -
python实现
# 导入模型包并创建实例
from docplex.mp.model import Model
model = Model(name='telephone_production')
# 创建待求解变量desk_production和cell_production
desk_production = model.continuous_var(name='desk_production')
cell_production = model.continuous_var(name='cell_production')
# 目标函数
model.maximize(12*desk_production+20*cell_production)
# 4个约束条件建议写约束条件都按照后两个的方式写等式
model.add_constraint(desk_production >= 100)
model.add_constraint(cell_production >= 100)
model.add_constraint(0.2*desk_production + 0.4*cell_production <= 400, ctname='ct_assembly')
model.add_constraint(0.5*desk_production + 0.4*cell_production <= 490, ctname='ct_painting')
# 打印模型信息
model.print_information()
# 对模型进行求解
sol = model.solve()
# 打印求解结果
print(sol)
- 输出结果
Model: telephone_production
- number of variables: 2
- binary=0, integer=0, continuous=2
- number of constraints: 4
- linear=4
- parameters: defaults
- objective: maximize
- problem type is: LP
solution for: telephone_production #这个显然是上面模型初始化里面给的名称telephone_production
objective: 20600
desk_production=300.000
cell_production=850.000
特殊情况——无解及解决方案
注意事项
官网上的例子用的是desk_production >= 1000这其实是不对的不方便后面处理无解的方法正确的应该是改cell_production好承接后文修改。且即使按照官网做的话本质上并没有修改一开始模型(m)的约束
c
t
−
a
s
s
e
m
b
l
y
=
m
.
a
d
d
−
c
o
n
s
t
r
a
i
n
t
(
0.2
∗
d
e
s
k
+
0.4
∗
c
e
l
l
<
=
400
)
ct_-assembly = m.add_-constraint( 0.2 * desk + 0.4 * cell <= 400)
ct−assembly=m.add−constraint(0.2∗desk+0.4∗cell<=400)并没有改成
c
t
−
a
s
s
e
m
b
l
y
=
m
.
a
d
d
−
c
o
n
s
t
r
a
i
n
t
(
0.2
∗
d
e
s
k
+
0.4
∗
c
e
l
l
<
=
1000
)
ct_-assembly = m.add_-constraint( 0.2 * desk + 0.4 * cell <= 1000)
ct−assembly=m.add−constraint(0.2∗desk+0.4∗cell<=1000)——因为无解的情况修改的是模型m的副本im。这一点根据官网给出修改后的方案结果可以验证。
即使直接在原始模型m上的约束条件改成
c
t
−
a
s
s
e
m
b
l
y
=
m
.
a
d
d
−
c
o
n
s
t
r
a
i
n
t
(
0.2
∗
d
e
s
k
+
0.4
∗
c
e
l
l
<
=
1000
)
ct_-assembly = m.add_-constraint( 0.2 * desk + 0.4 * cell <= 1000)
ct−assembly=m.add−constraint(0.2∗desk+0.4∗cell<=1000)也无法求得官方显示的正确结果。
因此本样例中选取的无解情况改的是cell_production。
无解的情况
显然把上述例子中cell_production值改大比如说cell_production >= 1000就得不到可行解此时结果返回为None
# 创建一个模型直接复制上面那个
im = model.copy()
# 因为模型是复制过来的已经没有上面那个desk_production变量了现在就是要把这个变量给拿出来重新换了个idesk的名字来使用
icell = im.get_var_by_name('cell_production')
# 加入新的约束——该约束实际上是覆盖了上面的那个desk_production >= 100的条件但是desk_production >= 100并没有删除
im.add_constraint(icell >= 1000)
# 输出一下此时的im模型信息
im.print_information()
# 如果没有得到可行解就返回None进入if去输出
ims = im.solve()
if ims is None:
print('- model is infeasible')
输出结果
对比上面的输出结果可以发现constraint的数量变多了说明原先的desk_production >= 100限制并没有删除所以说这种方法没办法删除constraint。而输出“- model is infeasible”表明该条件下确实没有可行解。
解决方案
- 一句话总结放宽限制约束硬约束转软约束
硬约束在任何情况都不能违反的约束上述的所有约束条件都是硬约束
软约束在某些情况可以违反的约束
- 如何转化
比如说想要把约束 0.2 d e s k − p r o d u c t i o n + 0.4 c e l l − p r o d u c t i o n ≤ 400 0.2desk_-production + 0.4cell_-production \le 400 0.2desk−production+0.4cell−production≤400 改为 0.2 d e s k − p r o d u c t i o n + 0.4 c e l l − p r o d u c t i o n ≤ 440 0.2desk_-production + 0.4cell_-production \le 440 0.2desk−production+0.4cell−production≤440
原始硬约束
0.2 d e s k − p r o d u c t i o n + 0.4 c e l l − p r o d u c t i o n ≤ 400 \\ 0.2desk_-production + 0.4cell_-production \le 400 0.2desk−production+0.4cell−production≤400
转化成软约束
0.2 d e s k − p r o d u c t i o n + 0.4 c e l l − p r o d u c t i o n ≤ 400 + o v e r t i m e \\ 0.2desk_-production + 0.4cell_-production \le 400 + overtime 0.2desk−production+0.4cell−production≤400+overtime
同时要限制overtime硬约束 o v e r t i m e ≤ 40 overtime \le 40 overtime≤40由于变成了软约束需要引入惩罚项。假定超过400后的额外费用为$2/h优化目标就需要修改为
m a x i m i z e : 12 ∗ d e s k − p r o d u c t i o n + 20 ∗ c e l l − p r o d c t i o n − 2 ∗ o v e r t i m e maximize: 12*desk_-production + 20*cell_-prodction - 2*overtime maximize:12∗desk−production+20∗cell−prodction−2∗overtime
对应代码直接在上述代码后面加
# 增加硬约束overtime
overtime = im.continuous_var(name='overtime', ub=40)
# 根据条件等式名称获得ct_assembly并对其等式右项进行修改
# 显然ct_assembly.rhs代表等式右项那么ct_assembly.lhs代表等式左项
ct_assembly = im.get_constraint_by_name('ct_assembly')
ct_assembly.rhs = 400 + overtime
# 增加惩罚项修改优化目标
im.maximize(12*desk_production + 20*cell_production - 2*overtime)
imsol = im.solve()
print(imsol)
输出结果
显然这是改了cell_production的取值范围与上文注意事项相互印证。
无界变量vs无界模型
- 定义与关系
无界变量变量没有上界||没有下界
无界模型模型的目标没有上界||没有下界说明限制条件少了
关系相互不影响
- Docplex中的默认设置
变量无上界下确界为0
如果不对上述例子中的3、4项约束解空间会是如下
解决线性规划问题LP的算法3种
单纯形法Simplex Optimizer
证明不细究主要是给一些结论证明可以参考这一篇具体怎么算可以参考这一篇
LP问题的解集是凸集
凸集的顶点一定是可行解
单纯形法核心思想找到每一个可行解带入目标函数后按照要求取最大/小值即可
几何解释由于解集是凸集只要找到顶点处也就是最极端的条件单纯形法本质是遍历凸集的每一个顶点
对偶单纯形法Dual-simplex Optimizer
对偶每个线性规划问题假设叫P问题都有一个关联的LP问题假设叫D问题那么P称之为原问题D称之为P的对偶
特点如果P是一个最小化问题则D是一个最大化问题反之亦然
对偶价格dual prices/影子价格对偶变量的值
怎么从原问题P推出对偶
原问题
m
i
n
i
m
i
z
e
c
T
x
s
u
b
j
e
c
t
t
o
A
x
≥
b
x
≥
0
minimize \quad c^Tx \\subject \quad to \quad \quad Ax \ge b \\x \ge 0
minimizecTxsubjecttoAx≥bx≥0
对偶
m
a
x
i
m
i
z
e
y
T
b
s
u
b
j
e
c
t
t
o
y
T
A
≤
c
T
y
≥
0
maximize \quad y^Tb \\subject \quad to \quad \quad y^TA \le c^T \\y \ge 0
maximizeyTbsubjecttoyTA≤cTy≥0
示例拿的这篇的
特性
- P问题中的约束都有一个关联的对偶变量 y i y_i yi
- D的任何可行解都是P的上限P的任何可行解都是D的下界
- LP问题中D和P的最终目标本质是等价的都是在边界处
- 如果原始问题很困难不便于解决可以考虑从对偶问题来解
障碍法Barrier Optimizer
证明暂时放过以后看了再补充
几何解释从可行区域的某个地方开始采用了预测-校正算法通过可行区域的中间路径不断调整路径。相比单纯形法沿着边缘走这个是直接从中间穿
调用方法
- 代码示例
# 修改lpmethod的值不同的值采用不同的算法lpmethod=4时采用的是用的障碍法求解
im.parameters.lpmethod = 4
# 求解并输出log
im.solve(log_output=True)
输出结果
- 参数取值
鉴于docplex是直接封装的cplex底层的值应该没改对应于docplex里面的lpmethod应该同cplex但是我在docplex里面没找到直接查找的参数说明要是有大佬知道麻烦踢我一脚不甚感激
>>> print(c.parameters.lpmethod.help())
method for linear optimization :
0 = automatic
1 = primal simplex
2 = dual simplex
3 = network simplex
4 = barrier
5 = sifting
6 = concurrent optimizers
参考链接
线性规划官方Docplex示例
Python 中 CPLEX 参数的帮助
| 阿里云国内75折 回扣 微信号：monov8 |
| 阿里云国际，腾讯云国际，低至75折。AWS 93折 免费开户实名账号 代冲值 优惠多多 微信号：monov8 飞机：@monov6 |