It is possible to deal with constraint optimization problem. Look at the following:
http://fenicsproject.org/documentation/dolfin/1.4.0/cpp/programmers-reference/nls/PETScSNESSolver.html
You set solver parameters and then define the nonlinear variational form of the problem along with appropriate lower and upper bounds in the function space that you are solving the problem. Finally, few lines can solve your problem (Assuming a_min and a_max are lower and upper bounds of the solution:
solver = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
a_min = Function(interpolate(Constant(0.0), S))
a_max = Function(interpolate(Constant(1.0), S))
(iter, converged) = solver.solve(a_min.vector(), a_max.vector())