cvxopt中关于不等式约束的错误（使用GLPK）

Question

我正在使用cvxopt来计算一个两人零和游戏的纳什均衡。

[-5, 3, 1, 8]
[ 5, 5, 4, 6]
[-4, 6, 0, 5]

这是我使用的代码（带有测试示例）。

from cvxopt import matrix, solvers
from cvxopt.modeling import op, dot, variable
import numpy as np


def solve_lp(a, b, c):
    """
    >>> a = matrix([[-5., 3., 1., 8., 1.],
    ...             [ 5., 5., 4., 6., 1.],
    ...             [-4., 6., 0., 5., 1.],
    ...             [-1.,-1.,-1.,-1., 0.],
    ...             [ 1., 1., 1., 1., 0.],
    ...             [-1., 0., 0., 0., 0.],
    ...             [ 0.,-1., 0., 0., 0.],
    ...             [ 0., 0.,-1., 0., 0.],
    ...             [ 0., 0., 0.,-1., 0.]])
    >>> b = matrix([0.,0.,0.,0.,1.])
    >>> c = matrix([0.,0.,0., 1.,-1.,0.,0.,0.,0.])
    >>> solve_lp(a, b, c)

    """
    variables = c.size[0]
    x = variable(variables, 'x')
    eq     =   (a*x == b)
    ineq   =   (x >= 0)
    lp = op(dot(c, x), [eq, ineq])
    lp.solve(solver='glpk')
    return (lp.objective.value(), x.value)

运行这段代码时出现了以下错误：

Traceback (most recent call last):
...
TypeError: 'G' must be a dense or sparse 'd' matrix with 9 columns

看起来cvxopt在处理ineq约束时抛出了一个异常，尽管我似乎遵循了建模示例中的约束语法。

我到目前为止尝试过的

通过将x乘以一个全是1的向量来修改代码：

def solve_lp(a, b, c):    
    variables = c.size[0]
    x = variable(variables, 'x')
    e = matrix(1.0, (1, variables))
    eq     =   (a*x == b)
    ineq   =   (e*x >= 0)
    lp = op(dot(c, x), [eq, ineq])
    lp.solve(solver='glpk')
    return (lp.objective.value(), x.value)

至少它能到达GLPK，但接着又出现了这个错误：

Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  8.000e+00  ratio =  8.000e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part = 6
*     0: obj =   0.000000000e+00  infeas =  0.000e+00 (0)
PROBLEM HAS UNBOUNDED SOLUTION
glp_simplex: unable to recover undefined or non-optimal solution

我该如何解决这个问题？

Answer 1

我觉得你应该参考这个网页上关于glpk求解器的用法：

https://github.com/benmoran/L1-Sudoku/blob/master/sudoku.py

如果你严格按照这个来做，你就能正确地修复并使用这个glpk求解器……

定义一个函数来解决问题：

def solve_plain_l1(A, b, solver='glpk'):

'''Find x with min l1 such that Ax=b,
using plain L1 minimization'''

n = A.size[1]

c0 = ones_v(2*n)

G1 = concathoriz(A,-A) # concatenate horizontally
G2 = concathoriz(-A,A)
G3 = -eye(2*n)
G = reduce(concatvert, [G1,G2,G3]) # concatenate vertically
hh = reduce(concatvert, [b, -b, zeros_v(2*n)])

u = cvxopt.solvers.lp(c0, G, hh, solver=solver)

v = u['x'][:n]

return v