0

我有一个相对简单的优化问题，但不是特别精通 scipy，无法弄清楚如何应用必要的约束。我的目标函数是最小化 10 元素向量y和x之间的绝对距离，其中x是 10x3 参数矩阵p的加权行和

x = p[:,0] + 2 * p[:,1] + 3 * p[:,2]

我需要将以下约束添加到参数矩阵p：

  1. p的每个元素，但是 >=0 和 <=1
  2. p的每一列之和必须为 1
  3. p的每一行之和不得超过1

我试图根据其他 SO 问题来实施约束，尽管我承认我并不完全理解它，也不完全理解它产生的错误。

# import packages
import numpy as np
from scipy.optimize import minimize

# create objective function
def objective(p, y):
    x = p[:,0] + 2*p[:,1] + 3*p[:,2]
    return np.sum(np.abs(y-x))

# define constraints
cons = [{'type':'ineq', 'fun': lambda x: x}, 
        {'type':'ineq', 'fun': lambda x: np.sum(x, 0)},   # row sum >= 0
        {'type':'ineq', 'fun': lambda x: 1 - np.sum(x, 0)},   # row sum <= 1
        {'type':'ineq', 'fun': lambda x: np.sum(x, 1) - 1},   # row sum >= 1
        {'type':'ineq', 'fun': lambda x: 1 - np.sum(x, 1)}]  # row sum <= 1

# generate target data (y)
y = np.array([2.48, 1.75, 1.13, 0.20, 0.20, 0.20, 0.0, 0.0, 0.0, 0.0])

# initialise starting param values
p0 = np.zeros((44,3))
p0[:,:] = 1/44

# solve
sol = minimize(objective, p0, method='SLSQP', constraints=cons)

运行此代码后，我收到以下错误消息：

AxisError：轴 1 超出维度 1 数组的范围

任何帮助深表感谢。

4

2 回答 2

2

这个答案基于 SuperKogito 的答案：

import numpy as np
from scipy.optimize import minimize

# create objective function
def objective(p, y):
    sep = len(y)
    x   = p[0:sep] + 2*p[sep:2*sep] + 3*p[2*sep:3*sep]
    return np.sum(np.abs(y-x))

def vec2mat(vec):
    sep = vec.size // 3
    return vec.reshape(sep, 3)

# define constraints
cons = [{'type':'eq', 'fun': lambda x: np.sum(vec2mat(x), axis=0) - 1},# sum cols == 1
        {'type':'ineq', 'fun': lambda x: 1-np.sum(vec2mat(x), axis=1)} # sum rows <= 1
        ]

# generate target data (y)
y = np.array([2.48, 1.75, 1.13, 0.20, 0.20, 0.20, 0.0, 0.0, 0.0, 0.0])

# initialise starting param values
p0      = np.zeros((10,3))
p0[:,0] = 0.001
p0[:,1] = 0.002
p0[:,2] = 0.003

# Bounds: 0 <= x_i <= 1
bnds = [(0, 1) for _ in range(np.size(p0))]

# solve
sol = minimize(objective, x0 = p0.flatten('F'), args=(y,), bounds=bnds, constraints=cons)
print(sol)

# reconstruct p
p_sol = vec2mat(sol.x)

给我

     fun: 1.02348743648716e-05
     jac: array([-1.,  1., -1.,  1., -1.,  1.,  1.,  1.,  1.,  1., -2.,  2., -2.,
        2., -2.,  2.,  2.,  2.,  2.,  2., -3.,  3., -3.,  3., -3.,  3.,
        3.,  3.,  3.,  3.])
 message: 'Optimization terminated successfully.'
    nfev: 1443
     nit: 43
    njev: 43
  status: 0
 success: True
       x: array([3.38285914e-01, 2.39989678e-01, 1.67911460e-01, 5.98105349e-02,
       4.84049819e-02, 1.44281667e-01, 1.42377791e-15, 1.74446343e-15,
       8.10053562e-16, 1.28295919e-15, 4.76418211e-01, 2.59584012e-01,
       2.20053270e-01, 5.22055787e-02, 2.00046857e-02, 1.43742204e-02,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       3.96292063e-01, 3.30281055e-01, 1.73992026e-01, 1.19261113e-02,
       3.71950067e-02, 8.98952507e-03, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00])

如果它满足所有约束，我们可以进一步检查我们的解决方案：

In [92]: np.all((p_sol >= 0) & (p_sol <= 1))
Out[92]: True

In [93]: np.sum(p_sol, axis=0)
Out[93]: array([1., 1., 1.])

In [94]: np.all(np.sum(p_sol, axis=1) <= 1)
Out[94]: True
于 2019-09-22T08:18:53.650 回答
1

正如@sascha提到的，scipy.minimize()只解决一维数组。因此，您需要展平矩阵（将其转换为向量），并在完成最小化后重建矩阵。这可以通过多种方式完成。在下面的代码中，我使用numpy.vstack()并numpy.concatenate()进行了转换。我还实施并加入了您的约束。请参阅以下代码：

# import packages
import numpy as np
from scipy.optimize import minimize

# create objective function
def objective(p, y, s):
    x   = p[0:sep] + 2*p[sep:2*sep] + 3*p[2*sep:3*sep]
    return np.sum(np.abs(y-x))

def vec2mat(vec):
    sep = int(vec.size / 3) 
    return np.vstack((np.vstack((vec[0:sep], vec[sep:2*sep])), vec[2*sep:3*sep]))

def mat2vec(mat):
    return np.concatenate((np.concatenate((mat[:,0], mat[:,1]), axis=0), mat[:,2]), axis=0)

def columns_sum(mat):
    return np.array([sum(x) for x in zip(*mat)])

def rows_sum(mat):    
    return np.array([sum(x) for x in mat])
def constraint_func(x):
    c1 = np.all(0 <= x) & np.all(x <= 1)         # 0 <= x_i <= 1
    c2 = np.all(rows_sum(vec2mat(x)) <= 1)       # row sum <= 1
    c3 = np.all(1 == columns_sum(vec2mat(x)))    # columns sum == 1
    if (c1 and c2 and c3) == True: return 1
    else                         : return 0

# define constraints
cons = [{'type':'eq', 'fun': lambda x: constraint_func(x) - 1 }]  

# generate target data (y)
y = np.array([2.48, 1.75, 1.13, 0.20, 0.20, 0.20, 0.0, 0.0, 0.0, 0.0])

# initialise starting param values
p0      = np.zeros((10,3))
p0[:,0] = 0.001
p0[:,1] = 0.002 
p0[:,2] = 0.003

# flatten p
p0_flat = mat2vec(p0)
sep     = int(p0_flat.size / 3) 

# solve
sol = minimize(objective, p0_flat, args=(y,sep), method='SLSQP', constraints=cons)
print(sol)

# reconstruct p
p_sol = vec2mat(sol.x)

不幸的是，尽管代码编译并正常工作，但它返回以下输出：

     fun: 5.932
     jac: array([-1., -1., -1., -1., -1., -1.,  1.,  1.,  1.,  1., -2., -2., -2.,
       -2., -2., -2.,  2.,  2.,  2.,  2., -3., -3., -3., -3., -3., -3.,
        3.,  3.,  3.,  3.])
 message: 'Singular matrix C in LSQ subproblem'
    nfev: 32
     nit: 1
    njev: 1
  status: 6
 success: False
       x: array([0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,
       0.001, 0.002, 0.002, 0.002, 0.002, 0.002, 0.002, 0.002, 0.002,
       0.002, 0.002, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003,
       0.003, 0.003, 0.003])

因此，您可能需要进一步调试代码（我怀疑约束部分），或尝试scipy.minimize_docs中提供的不同方法。

于 2019-09-19T17:02:46.060 回答