Vitaliy, that is interesting indeed. I wonder why the math for the periodic boundary conditions comes out as finite because one might think the power output of the light sources is constant, and if there is no transmission loss, then the amount of light in the room is as if the light has been running forever.

Here I tried adding reflections, although only for one reflection allowed, and a slider for the albedo.

Manipulate[ContourPlot[s@Total[Flatten[{f[x, y] @@@ pt, albedo {f[x - 1, y] @@@ pt, f[x + 1, y] @@@ pt, f[x, y - 1] @@@ pt, f[x, y + 1] @@@ pt}}]], {x, 0, 1}, {y, 0, 1}, PlotPoints -> 40, ColorFunction -> "IslandColors", Contours -> Table[min + c, {c, 0, .5, .1}]], {{pt, .5 + .3 Table[{Cos[2 Pi k/5], Sin[2 Pi k/5]}, {k, 5}]}, Locator, LocatorAutoCreate -> True}, {{s, Identity, "log scale"}, {Log -> "Yes", Identity -> "No"}}, {min, 13, 80}, {{albedo, .5}, 0, 1}, Initialization :> (f[x_, y_][a_, b_] := 1/((x - a)^2 + (y - b)^2);)]

It seems like this symmetric configuration is better with no albedo

pt = {{0.776`, 0.794`}, {0.20800000000000002`, 0.77`}, {0.226`, 0.226`}, {0.794`, 0.19`}, {0.53`, 0.512`}}

and this one is better with 50%

pt = {{0.773`, 0.864`}, {0.196`, 0.786`}, {0.226`, 0.226`}, {0.794`,0.15400000000000003`}, {0.646`, 0.516`}}

but I can't be sure, just from these pictures.