一个单一的线性规划是做不到的，因为你需要多个解决方案。
如果你的参与者变量少于（比如）30个，你可以做迭代LP，你把每个参与者的二进制赋值放在一个2幂的点积中，以产生一个规范顺序数字，逐步缩小规范顺序约束，以搜索更多的顺序答案。然而，对于高变量计数，这将遇到精度问题并停止工作。
这里是仍然使用迭代LP的解决方案。我最初猜测内部解搜索可以使用常用的丢番图求解器，但它不能，因为变量是有界的（二进制）。相反，我展示了一种方法，其中内部增量是二分搜索。请彻底测试这个，因为我没有。

from itertools import islice
from typing import Any, NamedTuple, Iterator

import numpy as np
from numpy.random import default_rng
from scipy.optimize import milp, LinearConstraint

class Problem(NamedTuple):
    skill_levels: np.ndarray
    salaries: np.ndarray
    max_salary: int = 20
    n_players: int = 2

    @classmethod
    def default_setup(cls) -> 'Problem':
        skill_levels = np.array([7, 13, 5, 1, 11])
        salaries = np.array([15,  8,  9, 12, 10])
        order = skill_levels.argsort()
        return cls(
            skill_levels=skill_levels[order],
            salaries=salaries[order],
        )

    @classmethod
    def big_setup(cls, m: int = 100_000) -> 'Problem':
        rand = default_rng(seed=0)
        skill_levels = rand.integers(low=1, high=20, size=m)
        salaries = rand.integers(low=1, high=20, size=m)
        order = skill_levels.argsort()
        return cls(
            skill_levels=skill_levels[order],
            salaries=salaries[order],
        )

    @property
    def n(self) -> int:
        return self.skill_levels.size

    @property
    def max_skill_level(self) -> float:
        return self.skill_levels.sum() + 0.5

def setup_milp(prob: Problem) -> tuple[
    dict[str, Any],  # milp kwargs
    LinearConstraint,  # skill sum constraint
]:
    skill_sum_constraint = LinearConstraint(A=prob.skill_levels, ub=prob.max_skill_level)

    return {
        # Maximize selected skill level
        'c': -prob.skill_levels,
        # Decision variables are integral
        'integrality': np.ones_like(prob.skill_levels, dtype=bool),
        # Decision variables are binary
        'bounds': (0, 1),

        'constraints': (
            # Exactly two players must be selected
            LinearConstraint(A=np.ones_like(prob.skill_levels), lb=prob.n_players, ub=prob.n_players),
            # Apply a fixed maximum salary
            LinearConstraint(A=prob.salaries, ub=prob.max_salary),
            # Evolving maximum for skill sum
            skill_sum_constraint,
        ),
    }, skill_sum_constraint

def outer_lp(milp_params: dict[str, Any], skill_sum_constraint: LinearConstraint) -> Iterator[int]:
    while True:
        res = milp(**milp_params)
        if not res.success:
            break
        skill_total = -round(res.fun)

        # res.x is one solution, but we ignore it and use the knowledge that skill_total produces at
        # least one solution.
        yield skill_total

        # ub is a read-only array, so replace it
        skill_sum_constraint.ub = np.array((skill_total - 0.5,), dtype=skill_sum_constraint.ub.dtype)

def log_traversal(prob: Problem, skill_total: int) -> Iterator[list[int]]:
    # for simplicity, hard-coded to two players
    first_limit = prob.n

    while True:
        first_idx = prob.skill_levels[:first_limit].searchsorted(skill_total)
        if first_idx >= first_limit:
            first_idx -= 1  # too high; use highest available
        first_limit = first_idx
        first_total = skill_total - prob.skill_levels[first_idx]

        second_idx = prob.skill_levels[:first_idx].searchsorted(first_total)
        if second_idx < first_idx:
            second_value = prob.skill_levels[second_idx]
            second_total = first_total - second_value
            if second_total == 0:
                # Handle duplicate values
                second_idx_arr, = (prob.skill_levels == second_value).nonzero()
                first_idx_arr = np.full_like(second_idx_arr, fill_value=first_idx)
                salary = prob.salaries[first_idx_arr] + prob.salaries[second_idx_arr]
                results = np.vstack((first_idx_arr, second_idx_arr))
                no_overlap = (first_idx_arr != second_idx_arr) & (salary <= prob.max_salary)
                salary = salary[no_overlap]
                results = results[:, no_overlap][:, salary.argsort()].T
                yield from results
        else:
            # if the second index needs to be higher than we can allow, there are no more solutions
            break

def all_solutions(prob: Problem) -> Iterator[tuple[
    int,  # skill total
    int,  # salary total
    list[int],  # player indices
]]:
    milp_params, skill_sum_constraint = setup_milp(prob)

    for skill_total in outer_lp(milp_params, skill_sum_constraint):
        for players in log_traversal(prob, skill_total):
            assert prob.skill_levels[players].sum() == skill_total
            yield skill_total, prob.salaries[players].sum(), players

def main() -> None:
    prob = Problem.big_setup()

    print('Skill, salary, assignments:')
    for skill, salary, players in islice(all_solutions(prob), 20):
        print(skill, salary, players)

if __name__ == '__main__':
    main()

字符串
与原始数据一起输出

Skill, salary, assignments:
24 18 [4 3]
18 17 [4 1]
16 19 [3 1]
14 20 [4 0]

型
使用this输出的代码中显示的较大数组来测试

Skill, salary, assignments:
38 12 [499999 497627]
38 12 [499999 475402]
38 12 [499999 492837]
38 12 [499999 484121]
38 12 [499999 486641]
38 12 [499999 479932]
38 12 [499999 494208]
38 12 [499999 486656]
38 12 [499999 496272]
38 12 [499999 494895]
38 12 [499999 486661]
38 12 [499999 475371]
38 12 [499999 484098]
38 12 [499999 489454]
38 12 [499999 475366]
38 12 [499999 475364]
38 12 [499999 475362]
38 12 [499999 496243]
38 12 [499999 480804]
38 12 [499999 486675]

型
n=20个记录。我没有非常小心的微优化，特别是围绕重复的东西，以便可以得到更快。

混音

对于1，000个玩家或更少的玩家来说，经历所有平分搜索对等麻烦没有多大意义。更改设计，使工资和技能矩阵保持为广播总和（对于旁边的1000个元素，每个仍然只占用4 MB的数量级）。并为LP添加第二个循环级别，以便它自然地优化技能第一，工资第二。这种方法要求技能和天赋都有独立的步骤。

from itertools import islice
from typing import Any, NamedTuple, Iterator

import numpy as np
from numpy.random import default_rng
from scipy.optimize import milp, LinearConstraint

class Problem(NamedTuple):
    skill_levels: np.ndarray
    salaries: np.ndarray
    skill_combos: np.ndarray   # for simplicity, hard-coded to two players
    salary_combos: np.ndarray
    max_salary: int = 20
    n_players: int = 2

    @classmethod
    def default_setup(cls) -> 'Problem':
        skill_levels = np.array([7, 13, 5, 1, 11])
        salaries = np.array([15,  8,  9, 12, 10])
        return cls.with_arrays(skill_levels, salaries)

    @classmethod
    def big_setup(cls, m: int = 1_000) -> 'Problem':
        rand = default_rng(seed=0)
        skill_levels = rand.integers(low=1, high=20, size=m)
        salaries = rand.integers(low=1, high=20, size=m)
        return cls.with_arrays(skill_levels, salaries)

    @classmethod
    def with_arrays(cls, skill: np.ndarray, salary: np.ndarray) -> 'Problem':
        return cls(
            skill_levels=skill,
            salaries=salary,
            skill_combos=np.tril(np.add.outer(skill, skill), k=-1),
            salary_combos=np.tril(np.add.outer(salary, salary), k=-1),
        )

    @property
    def n(self) -> int:
        return self.skill_levels.size

    @property
    def max_skill_level(self) -> float:
        return self.skill_levels.sum() + 0.5

def setup_milp(prob: Problem) -> tuple[
    dict[str, Any],  # milp kwargs
    LinearConstraint,  # skill sum constraint
    LinearConstraint,  # salary sum constraint
]:
    skill_sum_constraint = LinearConstraint(A=prob.skill_levels, ub=prob.max_skill_level)
    salary_sum_constraint = LinearConstraint(A=prob.salaries, lb=0)

    return {
        # Maximize selected skill level, and with less weight, minimize salary
        # Skill level must be integral for this to work
        'c': prob.salaries/2/prob.max_salary - prob.skill_levels,
        # Decision variables are integral
        'integrality': np.ones_like(prob.skill_levels, dtype=bool),
        # Decision variables are binary
        'bounds': (0, 1),

        'constraints': (
            # Exactly two players must be selected
            LinearConstraint(A=np.ones_like(prob.skill_levels), lb=prob.n_players, ub=prob.n_players),
            # Apply a fixed maximum salary
            LinearConstraint(A=prob.salaries, ub=prob.max_salary),
            # Evolving maximum for skill sum
            skill_sum_constraint,
            # Evolving minimum for salary sum
            salary_sum_constraint,
        ),
    }, skill_sum_constraint, salary_sum_constraint

def outer_lp(
    prob: Problem,
    milp_params: dict[str, Any],
    skill_sum_constraint: LinearConstraint,
    salary_sum_constraint: LinearConstraint,
) -> Iterator[tuple[int, int]]:
    while True:
        # Any salary
        salary_sum_constraint.lb = np.full_like(salary_sum_constraint.lb, 0)
        res = milp(**milp_params)
        if not res.success:
            break
        skill_total = round(prob.skill_levels.dot(res.x))
        salary_total = round(prob.salaries.dot(res.x))
        yield skill_total, salary_total

        # Fixed skill
        skill_sum_constraint.lb = np.full_like(skill_sum_constraint.lb, skill_total)
        skill_sum_constraint.ub = skill_sum_constraint.lb

        while True:
            # Salary floor
            salary_sum_constraint.lb = np.full_like(salary_sum_constraint.lb, salary_total + 0.1)
            res = milp(**milp_params)
            if not res.success:
                break

            salary_total = round(prob.salaries.dot(res.x))
            yield skill_total, salary_total

        # New skill ceiling
        skill_sum_constraint.lb = np.full_like(skill_sum_constraint.lb, 0)
        skill_sum_constraint.ub = np.full_like(skill_sum_constraint.ub, skill_total - 0.1)

def log_traversal(prob: Problem, skill: int, salary: int) -> np.ndarray:
    # for simplicity, hard-coded to two players
    yx = np.vstack((
        (prob.skill_combos == skill) &
        (prob.salary_combos == salary)
    ).nonzero())
    return yx.T

def all_solutions(prob: Problem) -> Iterator[tuple[
    int,  # skill total
    int,  # salary total
    list[int],  # player indices
]]:
    milp_params, skill_sum_constraint, salary_sum_constraint = setup_milp(prob)

    for skill, salary in outer_lp(prob, milp_params, skill_sum_constraint, salary_sum_constraint):
        for players in log_traversal(prob, skill, salary):
            assert prob.skill_levels[players].sum() == skill
            assert prob.salaries[players].sum() == salary
            yield skill, salary, players

def main() -> None:
    prob = Problem.big_setup()

    print('Skill, salary, assignments:')
    for skill, salary, players in islice(all_solutions(prob), 20):
        print(skill, salary, players)

if __name__ == '__main__':
    main()

Skill, salary, assignments:
38 3 [885 842]
38 4 [842 245]
38 4 [842 557]
38 4 [962 842]
38 5 [842 551]
38 5 [842 559]
38 5 [842 636]
38 5 [885 245]
38 5 [885 557]
38 5 [962 885]
38 6 [557 245]
38 6 [842 277]
38 6 [842 355]
38 6 [842 514]
38 6 [842 638]
38 6 [885 551]
38 6 [885 559]
38 6 [885 636]
38 6 [962 245]
38 6 [962 557]