History matching algorithms usually converge to the most prominent solution in the hypercube of parameter space defined by bound values. Here, we present a workflow to partition the parameter space into subdomains by defining a set of constraints. Then, a constrained history matching algorithm is developed to search each subdomain for a solution. This algorithm enables the engineers to solve the history matching problem subject to a set of general nonlinear/linear constraints on model parameters. The history matching problem definition follows a Bayesian framework, where the solution is obtained by maximizing the parameter's posterior probability density conditioned to the field data. With the proposed constrained algorithm, the optimization is subject to a set of constraints on model parameters. The optimizer is an iterative ensemble smoother and the constraints are enforced in a secondary update step at each optimization iteration by projecting the solutions to the feasible domain. The projection operator is derived from the Lagrangian form of the constrained problem, and based on linearizing the active set of constraints at the ensemble updates. The proposed constrained history matching algorithm and multi-solution search workflow are tested on an optimization test problem to validate its robustness and efficiency. Then history matching of a reservoir under water flooding is investigated where the history matching variables are the parameters for the relative permeability curves and the multipliers for the regional rock property fields. The constraints include relations between porosity and permeability multipliers as well as the relative permeability curve parameters. The constrained history matching algorithm could robustly find the feasible solutions which provided acceptable data matches. Moreover, with the application of the presented workflow, multiple solutions could be obtained for the history matching problem.