First of all, the important thing is that you need to consider at each iteration new cases where just one step to the right in one of the lists is made, e.g. starting from the solution $5+8+3$ you obtain $3+8+3\ ,\ 5+7+3\ ,\ 5+8+1$ . Never more than $1$ step and never in more than $1$ list otherwise you will get a lower sum.
So the problem is to deal with the "rollbacks" (steps to the left) in all the list except the one you go one step to the right. The thing so is to create an algorithm where you always compare all the possible rollbacks at each new choice of the current top sum.
You can do this incrementing by one each of the indices of the current top sum and adding this possible new top sums to the comparison.
I clarify with your example :
L1 = [5,3,3,1]
L2 = [8,7,0,0]
L3 = [3,1,1,1]
Start with indices 0,0,0 and best is 5+8+3 and so do not use this in the comparison,
use just the 3 incremented sums
(5+8+3,[0,0,0])
/ | \
(3+8+3,[1,0,0]) (5+7+3,[0,1,0]) (5+8+1,[0,0,1])
The current best now is 5+7+3 with indices 0,1,0 and do not consider that in the
next comparison, just use the other 5 sums
(3+8+3,[1,0,0]) (5+7+3,[0,1,0]) (5+8+1,[0,0,1])
/ | \
(5+0+3,[0,2,0]) (3+7+3,[1,1,0]) (5+0+3,[0,1,1])
The best is now 3+8+3 with indices 1,0,0 so we expand it ([0,0,1] has the same sum
so at the next iteration this sum will be taken)
As you see all the sums which compete to be the current best sum are taken into account at every new choice.
You can use a heap queue to implement the idea, remember to use a set where you add the new indices to compare so that to avoid duplicates in the queue. With $k$ lists, the time complexity is $O(Nk\cdot\max(\log k,\log N))$.
This is a Python implementation :
import heapq
def top_solutions(N,lists):
N_best = []
k, len_k_lists = len(lists), [len(x) for x in lists]
init_sol = (-sum(x[0] for x in lists),tuple(0 for x in range(k)))
heap_list = [init_sol]
seen = {init_sol[1]}
for _ in range(N) :
curr_best = heapq.heappop(heap_list)
N_best.append(curr_best)
ind = curr_best[1]
for x in range(k) :
if len_k_lists[x] > ind[x]+1 : r = tuple(c if i != x else c+1 for i,c in enumerate(ind))
else : continue
if r not in seen :
curr_sum = curr_best[0]-lists[x][r[x]]+lists[x][r[x]-1]
curr_value = (curr_sum,r)
heapq.heappush(heap_list,curr_value)
seen.add(r)
return [(-x[0],x[1]) for x in N_best]
N = 5
lists = [ [5,3,3,1],
[8,7,0,0],
[3,1,1,1] ]
print(top_solutions(N,lists))
output : [[16, [0, 0, 0]], [15, [0, 1, 0]], [14, [0, 0, 1]], [14, [0, 0, 2]], [14, [0, 0, 3]]]
Anyway you will obtain duplicate top sums, but with a different set of indices. If you are interested in just the values of the sums (e.g. you don't want two top sums equal to $14$ in the example ) then, first of all, you can eliminate the duplicate numbers in each list, and then check if you have $N$ different sums after $N$ iterations, if not you continue to run the procedure until you obtain them, I think there are also optimization to speed it up considerably. In this case, I think, the complexity is difficult to assess (anyway at worst roughly equal to the naive implementation, and for most inputs of $N$ less than it).