A variational approach to finite connectivity spin-glass-like models is developed and applied to describe the structure of optimal solutions in random satisfiability problems. Our variational scheme accurately reproduces the known replica symmetric results and also allows for the inclusion of replica symmetry breaking effects. For the 3-SAT problem, we find two transitions as the ratio α of logical clauses per Boolean variables increases. At the first one αs ≃ 3.96, a non-trivial organization of the solution space in geometrically separated clusters emerges. The multiplicity of these clusters as well as the typical distances between different solutions are calculated. At the second threshold αc ≃ 4.48, satisfying assignments disappear and a finite fraction B0 ≃ 0.13 of variables are overconstrained and take the same values in all optimal (though unsatisfying) assignments. These values have to be compared to αc ≃ 4.27, B0 ≃ 0.4 obtained from numerical experiments on small instances. Within the present variational approach, the SAT-UNSAT transition naturally appears as a mixture of a first and a second order transition. For the mixed 2+p-SAT with p < 2/5, the behavior is as expected much simpler: a unique smooth transition from SAT to UNSAT takes place at αc = 1/(1 − p).