You could solve this using Constraint Logic Programming (here ECLiPSe). You model the problem with variables ranging over a numeric domain, and invoke a built-in search routine. For simplicity, I have assumed length=weight.
:- lib(ic).
:- lib(branch_and_bound).
solve(Vars, Energy) :-
Vars = [TopSmall, TopLarge, BotSmall, BotLarge],
TopSmall :: 0..3, % how many small boxes on top
BotSmall :: 0..3, % how many small boxes on bottom
TopLarge :: 0..1, % how many large boxes on top
BotLarge :: 0..1, % how many large boxes on bottom
TopSmall + BotSmall #= 3, % total small boxes
TopLarge + BotLarge #= 1, % total large boxes
TopWeight #= TopSmall*1 + TopLarge*2, % total on top
BotWeight #= BotSmall*1 + BotLarge*2, % total on bottom
TopWeight #=< 3, % shelf capacities
BotWeight #=< 3,
Energy #= 2*TopWeight + 1*BotWeight, % Top shelf at double height
minimize(labeling(Vars), Energy). % find a minimal solution
% labeling(Vars). % alternatively, find all solutions