  We prove {the first} regularity theorem for the free boundary of solutions to
shape optimization problems involving integral functionals, for which the
energy of a domain $\Omega$ is obtained as the integral of a cost function
$j(u,x)$ depending on the solution $u$ of a certain PDE problem on $\Omega$.
The main feature of these functionals is that the minimality of a domain
$\Omega$ cannot be translated into a variational problem for a single (real or
vector valued) state function.
  In this paper we focus on the case of affine cost functions
$j(u,x)=-g(x)u+Q(x)$, where $u$ is the solution of the PDE $-\Delta u=f$ with
Dirichlet boundary conditions. We obtain the Lipschitz continuity and the
non-degeneracy of the optimal $u$ from the inwards/outwards optimality of
$\Omega$ and then we use the stability of $\Omega$ with respect to variations
with smooth vector fields in order to study the blow-up limits of the state
function $u$. By performing a triple consecutive blow-up, we prove the
existence of blow-up sequences converging to homogeneous stable solution of the
one-phase Bernoulli problem and according to the blow-up limits, we decompose
$\partial\Omega$ into a singular and a regular part. In order to estimate the
Hausdorff dimension of the singular set of $\partial\Omega$ we give a new
formulation of the notion of stability for the one-phase problem, which is
preserved under blow-up limits and allows to develop a dimension reduction
principle. Finally, by combining a higher order Boundary Harnack principle and
a viscosity approach, we prove $C^\infty$ regularity of the regular part of the
free boundary when the data are smooth.
