相当简单的LP构造;不过，您需要一个权重参数来确定平均误差相对于主要分配成本的重要性。尝试将下面的mean_err_weight在0.5和1.0之间更改，以观察差异。

import numpy as np
from scipy.optimize import milp, LinearConstraint
import scipy.sparse

def mean_assignment(
    asn_cost: np.ndarray,
    mean_cost: np.ndarray,
    target_mean: float,
    mean_err_weight: float = 1,
) -> tuple[
    np.ndarray,  # assignments
    float,  # mean
    float,  # error from mean
]:
    """
    Perform linear sum assignment with a secondary objective of minimizing deviation from a mean;
    in other words,
    minimize sum(asn_cost[assign]) + mean_err_weight*(sum(mean_cost[assign])/n - target_mean)

    :param asn_cost: Assigment costs, 2D
    :param mean_cost: Mean costs; must be same shape as assignment costs
    :param target_mean: Number to which the mean of mean_cost[assign] will be pulled
    :param mean_err_weight: Optimization objective weight for mean cost
    :return: Assignment array of same shape as asn_cost; resulting mean of assigned mean_cost; and
        resulting error from mean
    """

    if asn_cost.shape != mean_cost.shape:
        raise ValueError('Shape mismatch')

    # Ensure a long matrix during optimization
    transpose = asn_cost.shape[0] < asn_cost.shape[1]
    if transpose:
        asn_cost = asn_cost.T
        mean_cost = mean_cost.T

    m, n = asn_cost.shape

    # Variables:
    #    m*n assignments, binary
    #    mean deviance, continuous

    cost_coeff = np.concatenate((
        asn_cost.ravel(),    # Assignment costs
        (mean_err_weight,),  # Cost of error from mean
    ))
    integrality = np.ones(m*n + 1)  # Assignments are binary
    integrality[-1] = 0             # Mean error is continuous

    # There must be at most one assignment per row (Kronecker)
    row_excl = LinearConstraint(
        A=scipy.sparse.hstack(
            (
                scipy.sparse.kron(scipy.sparse.eye(m), np.ones((1, n))),
                scipy.sparse.csc_array((m, 1)),
            ),
            format='csc',
        ),
        lb=np.zeros(m),
        ub=np.ones(m),
    )

    # There must be exactly one assignment per column
    col_excl = LinearConstraint(
        A=scipy.sparse.hstack(
            (scipy.sparse.eye(n),)*m +
            (scipy.sparse.csc_array((n, 1)),),
            format='csc',
        ),
        lb=np.ones(n),
        ub=np.ones(n),
    )

    # The error from mean is taken as the absolute
    # err >= +sum(secondary[assign])/sum(assign) - target
    # err >= -sum(secondary[assign])/sum(assign) + target
    # cost_matrix*assign + ... - n*err <= n*target
    # cost_matrix*assign + ... + n*err >= n*target
    mean_err_abs = LinearConstraint(
        A=scipy.sparse.bmat([
            [mean_cost.ravel(), [-n]],
            [mean_cost.ravel(), [+n]],
        ], format='csc'),
        lb=[-np.inf, n*target_mean],
        ub=[n*target_mean, +np.inf],
    )

    result = milp(
        c=cost_coeff, integrality=integrality,
        constraints=(row_excl, col_excl, mean_err_abs),
    )
    if not result.success:
        raise ValueError(result.message)
    assignments = result.x[:-1].reshape((m, n)).astype(int)
    mean_err = result.x[-1]
    mean = mean_cost[assignments.astype(bool)].mean()
    if transpose:
        assignments = assignments.T
    return assignments, mean, mean_err

def test() -> None:
    assign, mean, mean_err = mean_assignment(
        asn_cost=np.array([
            [5, 2],
            [3, 4],
            [6, 7],
        ]),
        mean_cost=np.array([
            [1., 10.],
            [10., 1.],
            [2.,  2.],
        ]),
        target_mean=6.,
        mean_err_weight=1,  # change to 0.5 to see cost tradeoff
    )

    print(assign)
    print('mean =', mean)
    print('mean err =', mean_err)

if __name__ == '__main__':
    test()