Split the given array into K subsets such that maximum sum of all subsets is minimum

Use backtracking to check all the possible combinations and using a heuristic function you can reduce the search space to reach answer quickly. I've used one such solution to distribute the load as an array of n elements evenly across k partitions. Let me know if you need the code as well.

This is a pretty typical assignment problem that can be framed as a MILP optimisation. You need constraints for exclusive coverage, minimum cardinality, and sum-maximum correlation:

import time

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

a = [11, 11, 5, 10, 5]
n = len(a)
k = 2

# Variables: n assignments * k subsets, 1 maxsum

cost = np.concatenate((
    np.zeros(n*k), # assignments themselves do not have a cost
    np.ones(1),    # minimise subset maximum
))

integrality = np.concatenate((
    np.ones(n*k, dtype=np.uint8),        # assignments are binary
    np.zeros(1, dtype=np.uint8),   # maximum is continuous
))

bounds = Bounds(
    # Assignments are binary, maximum is unbounded
    lb=np.concatenate((np.zeros(n*k), np.full(1, fill_value=-np.inf))),
    ub=np.concatenate((np.ones(n*k), np.full(1, fill_value=+np.inf))),
)

# Each element must be assigned exactly once
coverage_constraint = LinearConstraint(
    A=sp.hstack((
        sp.kron(np.ones((1, k)), sp.eye_array(n)),  # assignment identity matrices
        sp.csc_array((n, 1)),  # max does not participate
    ), format='csc'),
    lb=1, ub=1,
)

# Each subset must be of cardinality at least one
cardinality_constraint = LinearConstraint(
    A=sp.hstack((
        sp.kron(sp.eye_array(k), np.ones((1, n))),  # assignment identity matrices
        sp.csc_array((k, 1)),  # max does not participate
    ), format='csc'),
    lb=1,
)

# The maximum must be at least the sum of each subset:
# a.I <= max   0 <= max - a.I
sum_constraint = LinearConstraint(
    A=sp.hstack((
        sp.kron(
            -sp.eye_array(k),  # assignment identity matrices
            np.reshape(a, (1, n)),  # array values to sum
        ),
        np.ones((k, 1)),
    ), format='csc'),
    lb=0,
)

t0 = time.perf_counter()
result = milp(
    c=cost, integrality=integrality, bounds=bounds,
    constraints=(
        coverage_constraint, cardinality_constraint, sum_constraint,
    ),
)
t1 = time.perf_counter()
assert result.success, result.message

print(result.message)
print(f'Time: {1e3*(t1-t0):.1f} ms')
assign, (max_sum,) = np.split(result.x, (-1,))
print(f'Maximum sum: {max_sum}')
assign = (assign > 0.5).reshape((k, n)).astype(np.uint8)
print('Assignments:')
print(assign)

Optimization terminated successfully. (HiGHS Status 7: Optimal)
Time: 3.6 ms
Maximum sum: 21.0
Assignments:
[[1 0 1 0 1]
 [0 1 0 1 0]]

This has the following sparsity pattern:

plt.imshow(
    sp.vstack((
        coverage_constraint.A, cardinality_constraint.A, sum_constraint.A
    )).toarray().clip(min=-1, max=+1),
)
plt.show()