问题描述

当我在下面的代码中运行test()时，我发现scipy.optimize.minimize()对于要求解的变量的不同初始化得到了不同的结果。

import numpy as np
from numpy.random import default_rng, Generator
from scipy.optimize import check_grad, minimize, Bounds, LinearConstraint, OptimizeResult

eta=10e-27
k=1e6#to balance the magnitude of the data

# The objective function
def sum_cost(
    f: np.ndarray,
    Adstr: float,
    coef1: np.ndarray,
    coef2: np.ndarray,
) -> float:
    return (k*coef1.dot(f**2) + coef2.dot(1/f) + Adstr)*1e5 #Scale results for easier optimization

# Minimize progress
def cvx1(
    A: np.ndarray,
    ds: np.ndarray,
    ts: np.ndarray,
    tr: np.ndarray,
    Cm: np.ndarray,
    eta: float = 10e-27,
) -> tuple[np.ndarray, OptimizeResult]:
    # parameter definition
    M, N = A.shape
    Adstr = (A/tr @ ds.T).sum()
    coef2 = (A*ts).ravel()
    coef1 = coef2 * eta

    # two constraints for my problem:
    # 1.x0>=0
    # 2.sum(x0,axis=1) <= Cm.reshape((M,)
    #  A      x         Cm
    # (M MN) (MN 1) <= (M 1)
    constraint = LinearConstraint(
        A=np.repeat(np.eye(M), N, axis=1),
        ub=Cm.ravel(),
    )

    # one way to initialize x0
    # Use analytic optimum when it obeys the constraint
    # Use equally divided Cm otherwise.
    optimum = (0.5 / (eta*k)) ** (1 / 3)
    print('optimum',optimum)
    x0 = np.zeros((M, N))
    for m in range(M):
        count = len(np.where(A[m, :] == 1)[0])
        # print(count)
        # print(C_m[0, m])
        x0[m, np.where((A[m, :]) == 1)[0]] = np.minimum(Cm[0, m] / count, optimum)
        x0[m, np.where((A[m, :] == 0))[0]] = 1e-16
    print('initial value for x0:',x0)

    result = minimize(
        fun=sum_cost,
        x0=x0.ravel(), # one way to initialize x0
        # x0=np.ones(M * N),# another way to initialize x0
        args=(Adstr, coef1, coef2),
        bounds=Bounds(lb=0),
        constraints=constraint,
        tol=1e-16,
        options={'maxiter': 1000},
    )
    f = result.x.reshape((M, N))

    print(result)

    if np.allclose(f, x0):
        print('Warning: f == initial; no optimization needed')

    return f, result

def test(M: int=3, N: int=6) -> None:
    rand = default_rng(seed=11)

    #some random parameters for the objective function

    # here just to produce a 0-1 matrix A:
    A = rand.uniform(low=0, high=1, size=(M, N))
    am_current = np.argmax(A, axis=0)
    a_current = np.zeros(A.shape)
    a_current[am_current, np.arange(A.shape[1])] = 1
    A=a_current
    print('A:', A)

    #other parameters
    ds = rand.uniform(500*8, 2_000*8, (1, N))
    ts = rand.uniform(10, 20, (1, N)) * ds
    tr = rand.uniform(50_000, 300_000, (M, N))
    C_m = rand.uniform(5e6, 20e6, (1, M))

    # solve
    x1, result1 = cvx1(A=A, ds=ds, ts=ts, tr=tr, Cm=C_m)

我运行了上面的代码两次，一次是使用如上所述定制的x0，另一次是使用x0所有1，但我得到了两个不同的结果，如下所示：

结果显示

1.如上所述定制的x0的结果

initial value for x0: [[1.0000000e-16 3.6840315e+06 1.0000000e-16 1.0000000e-16 1.0000000e-16
  3.6840315e+06]
 [1.0000000e-16 1.0000000e-16 3.6840315e+06 1.0000000e-16 1.0000000e-16
  1.0000000e-16]
 [3.6840315e+06 1.0000000e-16 1.0000000e-16 3.6840315e+06 3.6840315e+06
  1.0000000e-16]]
     fun: 86503.22680146241
     jac: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0.])
 message: 'Optimization terminated successfully'
    nfev: 19
     nit: 1
    njev: 1
  status: 0
 success: True
       x: array([1.0000000e-16, 3.6840315e+06, 1.0000000e-16, 1.0000000e-16,
       1.0000000e-16, 3.6840315e+06, 1.0000000e-16, 1.0000000e-16,
       3.6840315e+06, 1.0000000e-16, 1.0000000e-16, 1.0000000e-16,
       3.6840315e+06, 1.0000000e-16, 1.0000000e-16, 3.6840315e+06,
       3.6840315e+06, 1.0000000e-16])
Warning: f == initial; no optimization needed

============================== 1 passed in 0.18s ==============================

2. x0所有1的结果

fun: 87740.20159970068
     jac: array([ 0.        , -0.00097656,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        , -0.00097656,  0.        ,  0.        ,
       -0.00097656, -0.00097656,  0.        ])
 message: 'Optimization terminated successfully'
    nfev: 8336
     nit: 419
    njev: 417
  status: 0
 success: True
       x: array([2.57105896e+04, 3.42638038e+06, 1.51942392e+02, 1.51942383e+02,
       1.51942391e+02, 3.46166354e+06, 8.37974281e-01, 8.37972296e-01,
       4.25764167e+06, 8.37974124e-01, 8.37970035e-01, 8.37973655e-01,
       2.77064502e+06, 2.98881278e+06, 4.53859019e+04, 2.86306864e+06,
       3.19190058e+06, 3.59267809e+04])

========================= 1 passed, 1 warning in 0.43s =========================

一些疑问

1.如图所示，x0的两种不同初始化方式给予不同的结果。我的代码有什么问题吗？这些结果准确吗？
2.第一个定制的x0得到了更好的结果，不需要更多的优化（然而，它并不总是像我尝试的那样得到更好的结果），它已经是我的目标函数的一般最优解了吗？
3.顺便说一下，第一种初始化x0的方法被设计成遵循这样的规则，如果A[m][n]==0,x0[m][n]=0，但我不知道如何在Scipy中设置这种约束？
提前感谢和任何这些问题的答案是赞赏！