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)$$
j
(
u
,
x
)
=
-
g
(
x
)
u
+
Q
(
x
)
, where u is the solution of the PDE $$-\Delta u=f$$
-
Δ
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 $$
C
∞
regularity of the regular part of the free boundary when the data are smooth.