You can incorporate the area constraint with the Lagrange multiplier method; for example, the book of John Oprea Differential geometry and its applications has an accessible treatment of such variational problems.
However I recommend to look for a curve in parametric form $x=f(t), y=g(t)$ rather than in the graph form $y=f(x)$. The reason is as follows. The free part of the boundary (the curve you are looking for) wants to have constant curvature. Indeed, if a curve has different curvature $\kappa_1>\kappa_2$ at two points $p_1,p_2$, then we can push a little at the point $p_1$ to decrease the curvature, and push from the opposite side at $p_2$ to increase it. This perturbation will decrease the length without changing the area. A planar curve of constant curvature is a circular arc. So, the boundary wants to be a circular arc with given endpoints, bounding the given area. When such an arc is a graph $y=f(x)$, it is the solution. But in some cases (depending on your parameters) the extremal circular arc is not a graph, because it intersects some vertical line (either $x=x_1$ or $x=x_2$) twice. If you insist on the graph form of the solution, you create an obstacle problem, which makes the problem substantially more difficult.
On the geometric grounds (without doing the computations) I expect that when the obstacle comes into effect, there will not be a minimizer satisfying the given boundary conditions. Instead, a minimizing sequence will converge to a function which violates one or both boundary conditions, and whose graph is a circular arc tangent to the effective obstacle(s). For a simple example, take $x_1=-1, x_2=1$, prescribe the boundary conditions $f(-1)=f(1)=0$ and demand the area to be $10+\pi/2$. The minimizing functions will converge to $f(x)=5+\sqrt{1-x^2}$, which fails the boundary conditions. To recover the existence of minimizer, you can relax the problem by adding the length of vertical segments to the length of the proper graph: that is, the new variational problem is to minimize $|f(x_1)-y_1|+|f(x_2)-y_2|+\int_{x_1}^{x_2}\sqrt{1+f\,'(x)^2}\,dx$ subject to $\int_{x_1}^{x_2} f(x)\,dx=c$, and no boundary conditions. Still, an additional complication arises when the prescribed area $c$ is so small that $f$ wants to be negative.