The ray theory of light is equivalent to the Eikonal Equation, which in turn is essentially a slowly varying envelope approximation to Maxwell's equations. If we write the electric and magnetic field vectors as $\mathbf{E}\left(\mathbf{r}\right) = \mathbf{e}\left(\mathbf{r}\right) e^{i\,\varphi\left(\mathbf{r}\right)}$, $\mathbf{H}\left(\mathbf{r}\right) = \mathbf{h}\left(\mathbf{r}\right) e^{i\,\varphi\left(\mathbf{r}\right)}$ ($\mathbf{e}(\mathbf{r})$ and $\mathbf{h}(\mathbf{r})$ are the "envelope" functions) and assume a monochromatic field (so the implicit time dependence is $e^{-i\,\omega\,t}$, then the conditions for soundness of the Eikonal equations are:
$$\left|\mathbf{e}\right|^{-1} \left|\nabla \times\mathbf{e}\right| \ll \left|\nabla \varphi\right| \approx \left|k\right|$$
$$\left|\mathbf{h}\right|^{-1} \left|\nabla \times\mathbf{h}\right| \ll \left|\nabla \varphi\right| \approx \left|k\right|$$
You can very quickly tell that your $1{\rm \mu\,m}$ microlens is going to impose a phasefront on the field that wholesale violates the above conditions.
In short, ray otpics: forget it!
I show how to derive these conditions in my answers here and here.
You are in a lucky position, though, as your lenses are not much bigger than a wavelength. Therefore, full numerical simulations using finite element version of Maxwell's Equations are more than practicable, because you can mesh your system with elements well smaller than a wavelength and still have a simulation that will be swift. So here's how I would be going about the problem.
You MAY get some joy from simple analytical models grounded on Mie theory of scatterers. Born and Wolf, Chapter 14 has everything you need to know about Mie theory: your lenses, as nominal spheres are very much the kind of thing described by Mie theory of scattering. I would imagine a Mathematical workbook grounded on Mie theory as a rough guide and sanity check as you build more sophisticated simulations.
MEEP, an MIT baby (with a LISP interface - of course: what else?) for Maxwell equation solution by finite difference time domain simulation, is widespread and, once you get up to speed, would be a powerful softtool for your problem. I believe there is a Python interface for it, which would be much more fitting to the problem than LISP.