Multi-knapsack problem with aggregate objective function/objective with a soft limit

Working code follows, but here's how you would proceed with the formulation.

Formulation changes to the Multiple Knapsack problem given here

1. You will need two sets of new variables for each bin. Let's call them shortfall[j] (continuous) and filled[j] (boolean).

2. Shorfall[j] is simply the funding_requirement[j] - sum_i(funding[items i])

3. filled[j] is a Boolean, which we want to be 1 if the sum of the funding of each item in the bin is greater than its funding requirement, 0 otherwise.

4. We have to resort to a standard trick in Integer Programming that involves using a Big M. (A large number)

           if total_item_funding >= requirement, filled = 1
           if total_item_funding < requirement, filled = 0
   

This can be expressed in a linear constraint:

           shortfall + BigM * filled > 0

Note that if the shortfall goes negative, it forces the filled variable to become 1. If shortfall is positive, filled can stay 0. (We will enforce this using the objective function.)

5. In the objective function for a Maximization problem, you penalize the filled variable.

    Obj: Max sum(i,j) Xij * value_i + sum(j) filled_j * -100
   

So, this multiple knapsack formulation is incentivized to go close to each bin's funding requirement, but if it crosses the requirement, it pays a penalty.

You can play around with the objective function variables and penalities.

Formulation using Google-OR Tools

Working Python Code. For simplicity, I only made 3 bins.

from ortools.linear_solver import pywraplp


def create_data_model():
    """Create the data for the example."""
    data = {}
    weights = [48, 30, 42, 36, 36, 48, 42, 42, 36, 24, 30, 30, 42, 36, 36]
    values = [10, 30, 25, 50, 35, 30, 15, 40, 30, 35, 45, 10, 20, 30, 25]
    item_funding = [50, 17, 38, 45, 65, 60, 15, 30, 10, 25, 75, 30, 40, 40, 35]
    data['weights'] = weights
    data['values'] = values
    data['i_fund'] = item_funding
    data['items'] = list(range(len(weights)))
    data['num_items'] = len(weights)
    num_bins = 3
    data['bins'] = list(range(num_bins))
    data['bin_capacities'] = [100, 100, 80,]
    data['bin_funding'] = [100, 100, 50,]

    return data

def main():
    data = create_data_model()

    # Create the mip solver with the SCIP backend.
    solver = pywraplp.Solver.CreateSolver('SCIP')

    # Variables
    # x[i, j] = 1 if item i is packed in bin j.
    x , short, filled = {}, {}, {}
    for i in data['items']:
        for j in data['bins']:
            x[(i, j)] = solver.IntVar(0, 1, 'x_%i_%i' % (i, j))

    BIG_M, MAX_SHORT = 1e4, 500
    for j in data['bins']:
        short[j] = solver.NumVar(-MAX_SHORT, MAX_SHORT, 
                                 'bin_shortfall_%i' % (j))
        filled[j] = solver.IntVar(0,1,  'filled[%i]' % (i))

    # Constraints
    # Each item can be in at most one bin.
    for i in data['items']:
        solver.Add(sum(x[i, j] for j in data['bins']) <= 1)

    for j in data['bins']:
        # The amount packed in each bin cannot exceed its capacity.
        solver.Add(
            sum(x[(i, j)] * data['weights'][i]
                for i in data['items']) <= data['bin_capacities'][j])
        
        #define bin shortfalls as equality constraints
        solver.Add(
            data['bin_funding'][j] - sum(x[(i, j)] * data['i_fund'][i]
                for i in data['items']) == short[j])
        
        # If shortfall is negative, filled is forced to be true
        solver.Add(
            short[j] + BIG_M * filled[j] >= 0)


    # Objective
    objective = solver.Objective()

    for i in data['items']:
        for j in data['bins']:
            objective.SetCoefficient(x[(i, j)], data['values'][i])

    for j in data['bins']:
            # objective.SetCoefficient(short[j], 1)
            objective.SetCoefficient(filled[j], -100)

    objective.SetMaximization()

    print('Number of variables =', solver.NumVariables())
    status = solver.Solve()


    if status == pywraplp.Solver.OPTIMAL:
        print('OPTMAL SOLUTION FOUND\n\n')
        total_weight = 0
        for j in data['bins']:
            bin_weight = 0
            bin_value = 0
            bin_fund = 0

            print('Bin ', j, '\n')

            print(f"Funding {data['bin_funding'][j]} Shortfall \
            {short[j].solution_value()}")

            for i in data['items']:
                if x[i, j].solution_value() > 0:
                    print('Item', i, '- weight:', data['weights'][i], ' value:',
                          data['values'][i], data['i_fund'][i])
                    bin_weight += data['weights'][i]
                    bin_value += data['values'][i]
                    bin_fund += data['i_fund'][i]

            print('Packed bin weight:', bin_weight)
            print('Packed bin value:', bin_value)
            print('Packed bin Funding:', bin_fund)
            print()
            total_weight += bin_weight
        print('Total packed weight:', total_weight)
    else:
        print('The problem does not have an optimal solution.')


if __name__ == '__main__':
    main()

Hope that helps you move forward.