This is an optimization problem, where a function $f: \Bbb R^8 \rightarrow \Bbb R$ is defined as an area of the visibility polygon, corresponding to four points, situated inside a polygon with holes. So, your research should go in two directions:
Global optimization - a branch of numerical analysis, that deals with various iterative algorithms, able to find an argument of the optimized function, giving its minimum or maximum value over all the function domain. Sometimes the global optimization problem is solved as a sequence of local optimizations (with random starting points), however in general case there are no guarantee that the optimal solution will be found by this method. Please see this paper, where authors apply the Gradient descent method to maximize the area of the region, visible to the observer in a simple polygon (without holes).
Visibility polygon for a point $p$ in the plane configuration with possible boundary and obstacles is the polygonal region $V(p)$ of all points of the plane visible from the p. You'll need to compute visibility polygons many times for different points, so it makes sense to pre-process the original polygon (together with its holes) to make this computation faster. Please see this paper, where authors do exactly that - their approach needs $O(n^3log(n))$ pre-processing time, and it allows to compute the visibility polygon with $k$ vertices in $O((1+min(h,k))log(n)+k)$ time, where $h$ is the number of holes. The visibility polygon for four points will be a union of four visibility polygons (and this union can have holes as well) - so you'll need to compute this union. Please see Boolean operation on polygons as a starting point for your research as well.
Optimization problems on real numbers are intrinsically inexact and their solutions are approximations - the question is how close your solution should be to the optimal one. The gradient descent method might be enough in your case, where the bounding polygon is simply a rectangle, and holes also have simple shape.