Technological advances have led to the widespread use of computational models of increasing complexity, in both industry and everyday life. This helps to improve the design, analysis and operation of complex systems. Many computational models in the field of engineering consist of systems of coupled nonlinear partial differential equations (PDEs). As a result, optimization problems involving such models may lead to computational issues because of the large number of variables arising from the spatiotemporal discretization of the PDEs. In this work, we present a methodology for steady-state optimization, with nonlinear inequality constraints of complex large-scale systems, for which only an input/output steady-state simulator is available. The proposed method is efficient for dissipative systems and is based on model reduction. This framework employs a two-step projection scheme followed by three different approaches for handling the nonlinear inequality constraints. In the first approach, partial reduction is implemented on the equality constraints, while the inequality constraints remain the same. In the second approach an aggregation function is applied in order to reduce the number of inequality constraints and solve the augmented problem. The final method applies slack variables to replace the one aggregated inequality from the previous method with an equality constraint without affecting the eigenspectrum of the system. Only low-order Jacobian and Hessian matrices are employed in the proposed formulations, utilizing only the available black-box simulator. The advantages and disadvantages of each approach are illustrated through the optimization of a tubular reactor where an exothermic reaction takes place. It is found that the approach involving the aggregation function can efficiently handle inequality constraints while significantly reducing the dimensionality of the system.