找出 Pyomo 模型不可行的原因
Finding out reason of Pyomo model infeasibility
我得到了一个带有很多变量和约束的 pyomo 具体模型。
不知何故,我的模型中的一个变量违反了一个约束,这使得我的模型不可行:
WARNING: Loading a SolverResults object with a warning status into model=xxxx;
message from solver=Model was proven to be infeasible.
有没有办法询问求解器,不可行的原因?
因此,例如,假设我有一个名为 x
的变量,如果我定义以下 2 个约束,则模型将为 ofc。不可行。
const1:
x >= 10
const2:
x <= 5
以及我想要实现的目标,指出导致这种不可行性的约束和变量,以便我可以修复它。否则,对于我的大模型,很难找出导致这种不可行性的原因。
IN: write_some_comment
OUT: variable "x" cannot fulfill "const1" and "const2" at the same time.
许多求解器(包括 IPOPT)会在求解器终止时将变量值还给您,即使发现问题不可行。到那时,您确实有一些选择。
pyomo.util.infeasible
中提供的代码可能会对您有所帮助。 https://github.com/Pyomo/pyomo/blob/master/pyomo/util/infeasible.py
用法:
from pyomo.util.infeasible import log_infeasible_constraints
...
SolverFactory('your_solver').solve(model)
...
log_infeasible_constraints(model)
我不相信解算器在报告“不可行”后加载到模型中的任何数字。我认为任何求解器都不会保证这些数字的有效性。此外,除非包可以预测建模者的意图,否则不清楚它将如何列出不可行的约束。考虑 2 个约束条件:
C1: x <= 5
C2: x >= 10
X ∈ Reals, or Integers, ...
哪个是无效约束?这要看情况!重点是,根据求解器尝试的值来解开谜团似乎是一项不可能完成的任务。
一个可能的替代策略:用您认为有效的解决方案加载模型,并测试约束的松弛度。这个“加载的解决方案”甚至可以是一个空的情况,其中所有内容都归零(如果这在模型的上下文中有意义)。它也可以是一组通过单元测试代码尝试的已知可行解决方案。
如果您可以构建您认为有效的解决方案(忘记最优,只是一些有效的解决方案),您可以 (1) 加载这些值,(2) 遍历模型中的约束,(3)评估约束并寻找负松弛,以及 (4) 报告具有值和表达式的罪魁祸首
一个例子:
import pyomo.environ as pe
test_null_case = True
m = pe.ConcreteModel('sour constraints')
# SETS
m.T = pe.Set(initialize=['foo', 'bar'])
# VARS
m.X = pe.Var(m.T)
m.Y = pe.Var()
# OBJ
m.obj = pe.Objective(expr = sum(m.X[t] for t in m.T) + m.Y)
# Constraints
m.C1 = pe.Constraint(expr=sum(m.X[t] for t in m.T) <= 5)
m.C2 = pe.Constraint(expr=sum(m.X[t] for t in m.T) >= 10)
m.C3 = pe.Constraint(expr=m.Y >= 7)
m.C4 = pe.Constraint(expr=m.Y <= sum(m.X[t] for t in m.T))
if test_null_case:
# set values of all variables to a "known good" solution...
m.X.set_values({'foo':1, 'bar':3}) # index:value
m.Y.set_value(2) # scalar
for c in m.component_objects(ctype=pe.Constraint):
if c.slack() < 0: # constraint is not met
print(f'Constraint {c.name} is not satisfied')
c.display() # show the evaluation of c
c.pprint() # show the construction of c
print()
else:
pass
# instantiate solver & solve, etc...
报告:
Constraint C2 is not satisfied
C2 : Size=1
Key : Lower : Body : Upper
None : 10.0 : 4 : None
C2 : Size=1, Index=None, Active=True
Key : Lower : Body : Upper : Active
None : 10.0 : X[foo] + X[bar] : +Inf : True
Constraint C3 is not satisfied
C3 : Size=1
Key : Lower : Body : Upper
None : 7.0 : 2 : None
C3 : Size=1, Index=None, Active=True
Key : Lower : Body : Upper : Active
None : 7.0 : Y : +Inf : True
我得到了一个带有很多变量和约束的 pyomo 具体模型。
不知何故,我的模型中的一个变量违反了一个约束,这使得我的模型不可行:
WARNING: Loading a SolverResults object with a warning status into model=xxxx;
message from solver=Model was proven to be infeasible.
有没有办法询问求解器,不可行的原因?
因此,例如,假设我有一个名为 x
的变量,如果我定义以下 2 个约束,则模型将为 ofc。不可行。
const1:
x >= 10
const2:
x <= 5
以及我想要实现的目标,指出导致这种不可行性的约束和变量,以便我可以修复它。否则,对于我的大模型,很难找出导致这种不可行性的原因。
IN: write_some_comment
OUT: variable "x" cannot fulfill "const1" and "const2" at the same time.
许多求解器(包括 IPOPT)会在求解器终止时将变量值还给您,即使发现问题不可行。到那时,您确实有一些选择。
pyomo.util.infeasible
中提供的代码可能会对您有所帮助。 https://github.com/Pyomo/pyomo/blob/master/pyomo/util/infeasible.py
用法:
from pyomo.util.infeasible import log_infeasible_constraints
...
SolverFactory('your_solver').solve(model)
...
log_infeasible_constraints(model)
我不相信解算器在报告“不可行”后加载到模型中的任何数字。我认为任何求解器都不会保证这些数字的有效性。此外,除非包可以预测建模者的意图,否则不清楚它将如何列出不可行的约束。考虑 2 个约束条件:
C1: x <= 5
C2: x >= 10
X ∈ Reals, or Integers, ...
哪个是无效约束?这要看情况!重点是,根据求解器尝试的值来解开谜团似乎是一项不可能完成的任务。
一个可能的替代策略:用您认为有效的解决方案加载模型,并测试约束的松弛度。这个“加载的解决方案”甚至可以是一个空的情况,其中所有内容都归零(如果这在模型的上下文中有意义)。它也可以是一组通过单元测试代码尝试的已知可行解决方案。
如果您可以构建您认为有效的解决方案(忘记最优,只是一些有效的解决方案),您可以 (1) 加载这些值,(2) 遍历模型中的约束,(3)评估约束并寻找负松弛,以及 (4) 报告具有值和表达式的罪魁祸首
一个例子:
import pyomo.environ as pe
test_null_case = True
m = pe.ConcreteModel('sour constraints')
# SETS
m.T = pe.Set(initialize=['foo', 'bar'])
# VARS
m.X = pe.Var(m.T)
m.Y = pe.Var()
# OBJ
m.obj = pe.Objective(expr = sum(m.X[t] for t in m.T) + m.Y)
# Constraints
m.C1 = pe.Constraint(expr=sum(m.X[t] for t in m.T) <= 5)
m.C2 = pe.Constraint(expr=sum(m.X[t] for t in m.T) >= 10)
m.C3 = pe.Constraint(expr=m.Y >= 7)
m.C4 = pe.Constraint(expr=m.Y <= sum(m.X[t] for t in m.T))
if test_null_case:
# set values of all variables to a "known good" solution...
m.X.set_values({'foo':1, 'bar':3}) # index:value
m.Y.set_value(2) # scalar
for c in m.component_objects(ctype=pe.Constraint):
if c.slack() < 0: # constraint is not met
print(f'Constraint {c.name} is not satisfied')
c.display() # show the evaluation of c
c.pprint() # show the construction of c
print()
else:
pass
# instantiate solver & solve, etc...
报告:
Constraint C2 is not satisfied
C2 : Size=1
Key : Lower : Body : Upper
None : 10.0 : 4 : None
C2 : Size=1, Index=None, Active=True
Key : Lower : Body : Upper : Active
None : 10.0 : X[foo] + X[bar] : +Inf : True
Constraint C3 is not satisfied
C3 : Size=1
Key : Lower : Body : Upper
None : 7.0 : 2 : None
C3 : Size=1, Index=None, Active=True
Key : Lower : Body : Upper : Active
None : 7.0 : Y : +Inf : True