最大化二次 objective 函数服从 R 中的线性等式和不等式约束
Maximize quadratic objective function subject to linear equality and inequality constraints in R
我正在尝试寻找如何在 R 中同时使用等式和不等式约束来最大化二次函数：
最大化x' * H * x
受制于：Aeq * x = beq
A * x >= b
和 x >= 0
问题的简化版本可以是
最大化x^2 + y^2
受制于 x + y = 1
和 x, y >= 0
由于这是一个最大化问题，我无法使用 quadprog 包中的 solve.QP 函数。
我也尝试使用 constrOptim。但请注意，存在等式约束并且 constrOptim 需要在可行域内部进行初始猜测。因此 constrOptim 不能用于等式约束。
我也尝试在 alabama 包中使用 auglag。但我似乎没有得到最大化问题的正确答案。
如果问题是最小化问题，则简单问题的答案是 x = 0.5 和 y = 0.5。 auglag 和 solve.QP 都给了我这个答案。
但我正在寻找最大化问题的解决方案。几何学的答案是 (x = 1 and y = 0) OR (x = 0 and y = 1).
这不是一个完整的答案，但可能会向您展示一些替代算法方法。
问题
这个问题似乎是非凸的，这使得它很难解决（并且限制了可用的好软件的数量）。
一般非线性优化
正如 Christoph 在评论中提到的那样，一般非线性优化是一种可能的方法。当然，我们失去了关于全局最优解的保证。在内部使用优秀的开源软件 ipopt 将是一个很好的第一次尝试。
备选方案：Convex-Concave/Difference 凸面编程
您可能会考虑凸凹编程（它解决了一些全局性的简单问题），它有一个非常有效的启发式方法，称为 凸凹过程 (Yuille, Alan L., and Anand Rangarajan. "The concave-convex procedure." Neural computation 15.4 (2003): 915-936.)这应该比更一般的非线性方法更好。
我不确定在 R 中是否有很好的方法来做到这一点（无需手动完成），但在 Python 中有非常现代的开源研究库 dccp (based on cvxpy).
代码
from cvxpy import *
import dccp
from dccp.problem import is_dccp
x = Variable(1)
y = Variable(1)
constraints = [x >= 0, y >= 0, x+y == 1]
objective = Maximize(square(x) + square(y))
problem = Problem(objective, constraints)
print("problem is DCP:", problem.is_dcp())
print("problem is DCCP:", is_dccp(problem))
problem.solve(method='dccp')
print('solution (x,y): ', x.value, y.value)
输出
('problem is DCP:', False)
('problem is DCCP:', True)
iteration= 1 cost value = -2.22820497851 tau = 0.005
iteration= 2 cost value = 0.999999997451 tau = 0.006
iteration= 3 cost value = 0.999999997451 tau = 0.0072
('solution (x,y): ', 0.99999999872569856, 1.2743612156911721e-09)
Edit/Update
根据问题的大小（小），您还可以尝试 全局非线性求解器 ，例如 couenne。
代码
from pyomo.environ import *
model = ConcreteModel()
model.x = Var()
model.y = Var()
model.xpos = Constraint(expr = model.x >= 0)
model.ypos = Constraint(expr = model.y >= 0)
model.eq = Constraint(expr = model.x + model.y == 1)
model.obj = Objective(expr = model.x**2 + model.y**2, sense=maximize)
model.preprocess()
solver = 'couenne'
solver_io = 'nl'
stream_solver = True # True prints solver output to screen
keepfiles = True # True prints intermediate file names (.nl,.sol,...)
opt = SolverFactory(solver,solver_io=solver_io)
results = opt.solve(model, keepfiles=keepfiles, tee=stream_solver)
print("Print values for all variables")
for v in model.component_data_objects(Var):
print str(v), v.value
输出
Couenne 0.5.6 -- an Open-Source solver for Mixed Integer Nonlinear Optimization
Mailing list: couenne@list.coin-or.org
Instructions: http://www.coin-or.org/Couenne
Couenne: new cutoff value -1.0000000000e+00 (0.004 seconds)
NLP0012I
Num Status Obj It time Location
NLP0014I 1 OPT -0.5 6 0.004
Loaded instance "/tmp/tmpLwTNz1.pyomo.nl"
Constraints: 3
Variables: 2 (0 integer)
Auxiliaries: 3 (0 integer)
Coin0506I Presolve 11 (-1) rows, 4 (-1) columns and 23 (-2) elements
Clp0006I 0 Obj -0.9998 Primal inf 4.124795 (5) Dual inf 0.999999 (1)
Clp0006I 4 Obj -1
Clp0000I Optimal - objective value -1
Clp0032I Optimal objective -1 - 4 iterations time 0.002, Presolve 0.00
Clp0000I Optimal - objective value -1
NLP Heuristic: NLP0014I 2 OPT -1 5 0
no solution.
Clp0000I Optimal - objective value -1
Optimality Based BT: 0 improved bounds
Probing: 0 improved bounds
NLP Heuristic: no solution.
Cbc0013I At root node, 0 cuts changed objective from -1 to -1 in 1 passes
Cbc0014I Cut generator 0 (Couenne convexifier cuts) - 0 row cuts average 0.0 elements, 2 column cuts (2 active)
Cbc0004I Integer solution of -1 found after 0 iterations and 0 nodes (0.00 seconds)
Cbc0001I Search completed - best objective -1, took 0 iterations and 0 nodes (0.01 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
couenne: Optimal
"Finished"
Linearization cuts added at root node: 12
Linearization cuts added in total: 12 (separation time: 0s)
Total solve time: 0.008s (0.008s in branch-and-bound)
Lower bound: -1
Upper bound: -1 (gap: 0.00%)
Branch-and-bound nodes: 0
Print values for all variables
x 0.0
y 1.0
你可以用你可用的软件解决这个问题提供它可以处理非凸objectives。两点：
- 要绕过非线性等式约束的问题，您可以执行
x + y >= 1 和 x + y <= 1，这与 x + y = 1 的效果相同。这同样适用于一般 Aeq*x = b：Aeq*x <= b 和 Aeq*x >= b
评论中的- @jogo是对的：你只需要更改objective函数的符号即可。实际上，您正在 最大化 一个无界凸函数（即具有 [半] 正定 Hessian 矩阵），这仅在约束条件下才有意义。因此，同样适用于 最小化 无界凹函数（即具有 [半] 负定 Hessian 矩阵）。试想一下："going up as much as you can is the same as flipping the world upside down and going down as much as you can"
我正在尝试寻找如何在 R 中同时使用等式和不等式约束来最大化二次函数：
最大化x' * H * x
受制于：Aeq * x = beq
A * x >= b
和 x >= 0
问题的简化版本可以是
最大化x^2 + y^2
受制于 x + y = 1
和 x, y >= 0
由于这是一个最大化问题，我无法使用 quadprog 包中的 solve.QP 函数。
我也尝试使用 constrOptim。但请注意，存在等式约束并且 constrOptim 需要在可行域内部进行初始猜测。因此 constrOptim 不能用于等式约束。
我也尝试在 alabama 包中使用 auglag。但我似乎没有得到最大化问题的正确答案。
如果问题是最小化问题，则简单问题的答案是 x = 0.5 和 y = 0.5。 auglag 和 solve.QP 都给了我这个答案。
但我正在寻找最大化问题的解决方案。几何学的答案是 (x = 1 and y = 0) OR (x = 0 and y = 1).
这不是一个完整的答案，但可能会向您展示一些替代算法方法。
问题
这个问题似乎是非凸的，这使得它很难解决（并且限制了可用的好软件的数量）。
一般非线性优化
正如 Christoph 在评论中提到的那样，一般非线性优化是一种可能的方法。当然，我们失去了关于全局最优解的保证。在内部使用优秀的开源软件 ipopt 将是一个很好的第一次尝试。
备选方案：Convex-Concave/Difference 凸面编程
您可能会考虑凸凹编程（它解决了一些全局性的简单问题），它有一个非常有效的启发式方法，称为 凸凹过程 (Yuille, Alan L., and Anand Rangarajan. "The concave-convex procedure." Neural computation 15.4 (2003): 915-936.)这应该比更一般的非线性方法更好。
我不确定在 R 中是否有很好的方法来做到这一点（无需手动完成），但在 Python 中有非常现代的开源研究库 dccp (based on cvxpy).
代码
from cvxpy import *
import dccp
from dccp.problem import is_dccp
x = Variable(1)
y = Variable(1)
constraints = [x >= 0, y >= 0, x+y == 1]
objective = Maximize(square(x) + square(y))
problem = Problem(objective, constraints)
print("problem is DCP:", problem.is_dcp())
print("problem is DCCP:", is_dccp(problem))
problem.solve(method='dccp')
print('solution (x,y): ', x.value, y.value)
输出
('problem is DCP:', False)
('problem is DCCP:', True)
iteration= 1 cost value = -2.22820497851 tau = 0.005
iteration= 2 cost value = 0.999999997451 tau = 0.006
iteration= 3 cost value = 0.999999997451 tau = 0.0072
('solution (x,y): ', 0.99999999872569856, 1.2743612156911721e-09)
Edit/Update
根据问题的大小（小），您还可以尝试 全局非线性求解器 ，例如 couenne。
代码
from pyomo.environ import *
model = ConcreteModel()
model.x = Var()
model.y = Var()
model.xpos = Constraint(expr = model.x >= 0)
model.ypos = Constraint(expr = model.y >= 0)
model.eq = Constraint(expr = model.x + model.y == 1)
model.obj = Objective(expr = model.x**2 + model.y**2, sense=maximize)
model.preprocess()
solver = 'couenne'
solver_io = 'nl'
stream_solver = True # True prints solver output to screen
keepfiles = True # True prints intermediate file names (.nl,.sol,...)
opt = SolverFactory(solver,solver_io=solver_io)
results = opt.solve(model, keepfiles=keepfiles, tee=stream_solver)
print("Print values for all variables")
for v in model.component_data_objects(Var):
print str(v), v.value
输出
Couenne 0.5.6 -- an Open-Source solver for Mixed Integer Nonlinear Optimization
Mailing list: couenne@list.coin-or.org
Instructions: http://www.coin-or.org/Couenne
Couenne: new cutoff value -1.0000000000e+00 (0.004 seconds)
NLP0012I
Num Status Obj It time Location
NLP0014I 1 OPT -0.5 6 0.004
Loaded instance "/tmp/tmpLwTNz1.pyomo.nl"
Constraints: 3
Variables: 2 (0 integer)
Auxiliaries: 3 (0 integer)
Coin0506I Presolve 11 (-1) rows, 4 (-1) columns and 23 (-2) elements
Clp0006I 0 Obj -0.9998 Primal inf 4.124795 (5) Dual inf 0.999999 (1)
Clp0006I 4 Obj -1
Clp0000I Optimal - objective value -1
Clp0032I Optimal objective -1 - 4 iterations time 0.002, Presolve 0.00
Clp0000I Optimal - objective value -1
NLP Heuristic: NLP0014I 2 OPT -1 5 0
no solution.
Clp0000I Optimal - objective value -1
Optimality Based BT: 0 improved bounds
Probing: 0 improved bounds
NLP Heuristic: no solution.
Cbc0013I At root node, 0 cuts changed objective from -1 to -1 in 1 passes
Cbc0014I Cut generator 0 (Couenne convexifier cuts) - 0 row cuts average 0.0 elements, 2 column cuts (2 active)
Cbc0004I Integer solution of -1 found after 0 iterations and 0 nodes (0.00 seconds)
Cbc0001I Search completed - best objective -1, took 0 iterations and 0 nodes (0.01 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
couenne: Optimal
"Finished"
Linearization cuts added at root node: 12
Linearization cuts added in total: 12 (separation time: 0s)
Total solve time: 0.008s (0.008s in branch-and-bound)
Lower bound: -1
Upper bound: -1 (gap: 0.00%)
Branch-and-bound nodes: 0
Print values for all variables
x 0.0
y 1.0
你可以用你可用的软件解决这个问题提供它可以处理非凸objectives。两点：
- 要绕过非线性等式约束的问题，您可以执行
x + y >= 1和x + y <= 1，这与x + y = 1的效果相同。这同样适用于一般Aeq*x = b：Aeq*x <= b和Aeq*x >= b
评论中的 - @jogo是对的：你只需要更改objective函数的符号即可。实际上，您正在 最大化 一个无界凸函数（即具有 [半] 正定 Hessian 矩阵），这仅在约束条件下才有意义。因此，同样适用于 最小化 无界凹函数（即具有 [半] 负定 Hessian 矩阵）。试想一下："going up as much as you can is the same as flipping the world upside down and going down as much as you can"