Model order reduction with Galerkin projection applied to nonlinear optimization with infeasible primal‐dual interior point method

It is not new that model order reduction (MOR) methods are employed in almost all fields of engineering to reduce the processing time of complex computational simulations. At the same time, interior point methods (IPMs), a method to deal with inequality constraint problems (which is little explored in engineering), can be applied in many fields such as plasticity theory, contact mechanics, micromechanics, and topology optimization. In this work, a MOR based in Galerkin projection is coupled with the infeasible primal‐dual IPM. Such research concentrates on how to develop a Galerkin projection in one field with the interior point method; the combination of both methods, coupled with Schur complement, permits to solve this MOR similar to problems without constraints, leading to new approaches to adaptive strategies. Moreover, this research develops an analysis of error from the Galerkin projection related to the primal and dual variables. Finally, this work also suggests an adaptive strategy to alternate the Galerkin projection operator, between primal and dual variable, according to the error during the processing of a problem.

of the 20th century with the development of the infeasible primal-dual IPM, which is explained in other works [17][18][19][20][21][22][23][24] ; such variation ensures much more stability for the simulation. In some fields of computational simulations in engineering, such as topology optimization 25,26 and mechanics, 27,28 the infeasible primal-dual IPM already has applications; however, in the opinion of the authors, the method is still far from its true potential.
Because of that, the proposal of this work is a MOR developed from the Galerkin projection and the infeasible primal-dual IPM. The originality of this work is the application of the infeasible primal-dual IPM with special variables, with the Galerkin projection in a particular way; such configuration permits the application of a MOR with the processing similar to problems without constraints, leading to new approaches to adaptive strategies. Besides the originality of this work, the authors are concerned to demonstrate the potential and limitations of such combination of methods; and this way, permitting the applicability in fields like plasticity theory, contact mechanics, multiscale, and topology optimization to future research works, which can use this work as theoretical support and orientation. Based on all explained so far, and now using a more detailed approach, this work has the following objectives.
1. Develop a MOR to a multiconstraint problem with inequality. This MOR uses the Galerkin projection applied in only one field per time of the problem; then, because of the formulation of the infeasible primal-dual IPM, a Schur complement of this field is applied. The combination of these methods results in a MOR that permits to deal with an adaptive strategy that works equivalent in the processing as a problem without constraint. 2. Analysis of error from the primal and dual field. With the projection isolated in one field per time, the error of the respective Galerkin projection is analyzed, discriminating each parcel related to the primal or dual field. As expected, this work gives a special attention related to the parcel associated with the Karush-Kuhn-Tucker (KKT) condition. 3. The development of an adaptive strategy based on previous objectives. Now that the projected field is isolated and its error is analyzed, it is possible to propose an adaptive strategy based on the selection of the reduced field and the respective Schur complement, which combination minimizes the processing time. Two examples are tested with variations of the selection of the reduced field and the Schur complement; these examples not only have the objective to demonstrate the time reducing but also show the dependence between the error and the processing time.
After defining the objectives, the development of sections follow them. The second section introduces the idea behind the IPM; it explained the development of the infeasible primal-dual IPM, from the concept of the barrier function, passing through the classical infeasible primal-dual IPM, to the general approach with special variables, as described in the work of Svanberg. 26 After the method is formulated, it is applied to the Schur complement related to the primal and dual variable. This section ends with an explanation about the IPM algorithm and the recommended adaptations of the method to its application.
The third section explains the first objective of this research; until now, the second section was only a preparation to understand this section. The third section described the general idea of the Galerkin projection. Then, it develops the formulation related to the Galerkin projection in one field per time and its respective Schur complement to the primal and dual field. Like the previous one, this section ends with an explanation about the IPM algorithm with MOR and the recommended adaptations to its application.
The fourth section describes the second objective of this research; the Galerkin error is analyzed with the formulations obtained from the second and third sections. In these Galerkin errors, from the primal and dual variable, the parcels of each error are analyzed individually. Sequentially, a Galerkin projection with both fields simultaneously is demonstrated, and a combination of both errors and the dependency between them are analyzed. Moreover, the last part of a sensibility analysis is developed.
The fifth and sixth sections, together, conclude the third objective. The fifth section describes the adaptive strategy in a complete online analysis, with many applications extracted from the works of Nigro et al, 1,2 but, now, the MOR enters like the accelerator of the preconditioned conjugate gradient (PCG). This section starts from the definition of the selection of the snapshot and ends to the adaptive strategy to the selection of the Galerkin projection. Finally, these strategies are applied to the problems described in the sixth section. In the sixth section, the numerical examples are tested with different selections of Galerkin projections as well as different Schur complements. In this way, the processing time and Galerkin errors measurements provide relevant conclusions related to the advantages and disadvantages of the MOR develop in this work as well as possible applications for numerical simulations in engineering.
Before starting the next sections, it is important to explain that, to resume the formulation, this work uses some definitions of vector and matrices, ie,

Classical infeasible primal-dual IPM
The MOR described in this work starts with the definition of the optimization problem. It has chosen a nonlinear optimization with multiple inequality-constrained and primal variable, or design variable "x", One approach to deal with inequality constraints is to apply a Lagrangian function and to use a Newton-Raphson method. In inequality constraint, the problem with this approach has an additional algorithm to detect active constraints; such algorithm changes constantly the dimension of the system and generates a considerable instability and loss of convergence. This additional algorithm has a very common application in engineering simulation, especially in a problem that involves contact mechanics and plasticity theory. Another approach, which is applied in this work, is the addition of slack variables to become from inequalities to equalities, ie, The slack variables named "s" must have positive values to ensure that they fulfill the remaining value to the equality. To deal with this property of the slack variable, it is applied a continuous, differentiable, and smooth approach to ensure a limitation in the domain. Such approach is called barrier function. The barrier function is a continuous function that has the behavior to reach a high value when close to the limit of the feasible region. Based on these modifications, the optimization problem becomes The natural logarithm function described on (3) is the barrier function. It is possible to realize that, when the value inside the logarithm function reaches a value close to the defined limit (in this case, zero), it results in a high negative value, ie, However, there is the barrier parameter " " to ensure, during the simulation, that system (3) implies to system (2), ie, Another aspect to consider related to the slack variables and barrier function is that they belong to the infeasible version of the IPM, as applied in related works, [18][19][20][21] which permits the violation, eventually, of the inequality constraint. To solve the optimization problem described in (3), it develops a Lagrangian function In the Lagrangian function, it is possible to realize that the slack variable can replace the value of the constraint function to satisfy the KKT condition, in which the value is imposed to be equal to zero. Such imposition is obtained, in steps, by the decrement of the barrier parameter The barrier parameter decreases during the simulation of the problem, ensuring the respect of the KKT condition at the end of the resolution.

Infeasible primal-dual IPM with artificial variables
In many cases, to support the application of slack variables and barrier function, the optimization problem as described in (1) has some instabilities and lack of generality. To handle with these limitations, the modification described in the work of Svanberg 26 is implemented, ie, The problem described in (8) is similar to (1), but it has advantages caused by artificial variables named "y" and "z". The artificial variables intend to ensure stability and generality of the problem. The artificial variable "y" has the property to ensure at least one point that satisfies the KKT condition; the artificial variable "z" permits a generalization of the optimization problem, which includes minimizing the maximum (min-max) and least square problems. This generalization of the optimization problem is very useful in some engineering simulation fields, such as topology optimization. Another aspect to consider is that the problem described in (8) is applied with specific convex approximation functions described in the works of Svanberg 25,26 ; however, the formulation is so general that it can be applied with different locally convex functions.
Still analyzing the definition of the optimization problem, it is possible to notice that, to satisfy the KKT condition, the artificial variables must have values higher than zero. In addition, in this work and in most cases, the optimization must be made in a limited domain to respect parameter of design; to ensure that, the design variable "x" has values between " " and " ". To solve these domain limitation problems, there are added functions in the optimization problem (8) to penalize these variables when they are close to the limit of the domain, ie, Problem (9) is an extension of problem (3); the idea is to ensure more stability and increase the range of application to the classical problem. Now, to start the interior method development, it is necessary to formulate a Lagrangian function from (9), ie, Now that the Lagrangian function is defined, it is necessary to find the critical point or the solution to the problem. The solution is developed by the Newton-Raphson method using primal and dual variables. Before starting the Newton-Raphson method, it is necessary to linearize the space. The linearization of the space results that each iteration of the optimization problem is developed at one point of the manifold , where the manifold  is defined by In this manifold, a point is a tangent vector space  u , and, consequently, The geometrical representation is shown in Figure 1. Figure (1) shows the vector space  u  amplified to a prism. It is important to notice that the representation by a prism is unusual (vector spaces are represented by planes in the manifold) but the additional dimension plays a useful role in the following sections of this work. Now, in the vector space  u , the fields are linearized, ie, In the same vector space, developing Taylor's series expansions from the Lagrangian function (10), the terms of second and higher order, as the first order null terms, are ignored. As a result, Taylor's series expansions in the critical point, arranged in the matrix form, are

FIGURE 1
The vector space of the manifold To simplify the formulation, it is defined as It is possible to realize in (15) that the fourth row (top to bottom) of (14) is multiplied by the slack operator defined in (7); the objective is the simplification of the system further. The gradients from (15) are In (16), the matrices shown previously are Moreover, the Hessians (14) are As noticed, the system of (14) is not symmetric, but with the manipulation of the last two rows, it is possible to reduce the dimension and to ensure symmetry. Isolating the increment of the slack variable in the fourth row, and replacing in the fifth row, the new system is where In the second row of system (25), it is possible to obtain directly the increment of the special variable "y", ie, Moreover, the increment of the slack variable is obtained directly of the fourth row from (14), ie, Matrix (25) can be considered as the curvature of a saddle function in this case because all functions related to each domain are forced to be convex. An approach that allowed to treat saddle function as convex function, at least locally, is the Schur complement. The Schur complement reduces the dimension of the system to the dimension of one domain; the price to pay is the calculation of one inverse matrix. However, such price is mitigated if the inverted matrix is a diagonal matrix. Supported by this idea, this work develops the Schur complement in two domains: one related to the primal variable and another related to the dual variable.

Schur complement related to the primal variable
This approach of the Schur complement is indicated when the number of constraints (dual variables) is greater than the number of design variables (primal variables). It is also indicated when the inverted Hessian of the design variable is expensive. To develop this Schur complement, the increment of the artificial variable "y" is necessary as described in (27), and replacing in the third row of the system (25), ie, The new arrangement of the system (25) becomes Based on this arrangement and using the Schur complement as detailed in Appendix A, results in the system where The new system results Furthermore, the increment of the dual variable is found from the second row of system (30), ie,

Schur complement related to the dual variable
This approach of the Schur complement is indicated when the number of design variables (primal variables) is greater than the number of constraints (dual variables). It is not indicated when the inverted Hessian of the design variable is expensive. To develop this Schur complement, it is necessary to separate system (25) in a specific subsystem of the three first rows, ie, Based in this subsystem, using the Schur complement as detailed in Appendix A and regrouping with the last row of the artificial variable "z" on system (25) It is possible to find The new system results In addition, the increment of the primal variable is found from the first row of system (25), ie,

The IPM algorithm and KKT conditions
Heretofore, the IPM algorithm related to the described problem is explained. As described in the last two subsections, the procedures to update all the variables using the Schur complement can be resumed in Algorithm 1A and Algorithm 1B from system (25).
When all variables are updated, the KKT condition must be checked. In this case, it is possible to check this condition by the residuals. In the analysis, the residual is calculated from the gradients (16), (17), (18), (19), and (20), ie, When the residual is defined, the IPM algorithm can be developed. In the algorithm applied in this work, the barrier parameter decreases every time that the norm of residuals described in (48) reaches values below the value of the current barrier parameter. Lastly, when the barrier parameter reaches the minimum specified value, it is considered that the KKT condition is satisfied, resulting in the solution. For a better understanding, a scheme of the IPM resolution is shown in Flowchart 1 in the next section.

Global Galerkin projection
In an optimization with constraint, it is necessary to handle with at least two fields, represented by the design variable (primal variable) and the Lagrange multiplier (dual variable). In a common MOR, all the fields are reduced, as described in the works of Xiao et al 9 and Zahr and Farhat. 13 Another approach, as described in the works of Kaulmann et al, 5,6 is a localized MOR, but an offline analysis previously is necessary. In some cases, the complete online MOR with several fields can reach a similar processing time when compared to the full model.
To start an online MOR from infeasible primal-dual IPM, the proposal is to work with the reduction, through the Galerkin projection, in only one field per time. To develop that, the application of the Schur complement to achieve this goal is shown. In the type of problem described in the previous section, this is possible with low numerical effort. During the processing of the problem, the update and the error control of the reduced model can be treated as an optimization without constraint, as described in the work of Nigro et al. 1,2 The application of the Schur complement, in the point of view of the size of the domain, can be considered as a reduction of dimension. The Schur complement reduces the dimension of the system to the dimension of one field (usually the smaller) during the solution. In a model with two fields, when dimensions are in the same order of magnitude, the selection of the Schur complement is not so effective. Now, presenting the infeasible primal-dual IPM and the Schur complement, aims in the MOR are necessary. The objective of a MOR, based in the Galerkin projection, is to find a vector space  u such that The vector space  u is the reduced space, and it can be represented by the plane inside the prism, ie, the vector space  u . To develop this MOR based in the Galerkin projection, it is necessary to start with the matrix and vectors from (25). In the Galerkin projection, the increment of the domain is decomposed in two additive terms, ie, the range and null space, Vector Δu Ψ represents the range space operator; it always belongs to the space vector  u . Vector Δu Ω represents the null space; it represents the error of Δu Ψ related to Δũ, and it is always orthogonal to the space vector  u . Figure 2 shows the reduction of dimension, and the reason to the point of the manifold  be represented by a prism, as declared in the previous section. Back to the explanation of (50), matrix represents the global Galerkin operator, and matrix K represents the global null space operator. Each global operator belongs to complementary spaces To calculate the error associated to the projection, the system described in (25) is operated by T in both sides, and considering the decomposition of (50), The first term on the left side in Equation (52) represents the part of the gradient that belongs to the range space and provides the result of the MOR; the second term on the left side in Equation (52) represents the part of the gradient that belongs to the null space, and, consequently, does not belong to the MOR; this last term contains the error of the reduced space related to the full space.
Related to (25), the global Galerkin operator can obtain the projection until three domains. However, it is possible to notice that the simultaneous projections generate cross terms; the cross terms difficult a more accurate analysis of the error related to each domain. Because of that, one of the proposals of this work is the application of the projection in one domain per time (in addition to the error being measured more precisely) and there is an advantage related to the operation with the Schur complement that is demonstrated in the next two sections.

MOR of the primal field
In this section, the main hypothesis (a reduction of one field) considers the decomposition of the design variable (primal variable), ie, where x is the local range space operator related to the primal space, and x is the local null space operator. The main variables are arranged in the following decomposition: When this global Galerkin operator and the global null space operator are replaced in (52), the parcel of the error is represented by the term Term (56) is useful to measure the error related to the full model, including the application in an adaptive strategy to reduce the error of the MOR. Finally, the terms related to the projection are Applying the Schur complement as developing in (31) and (37), where As demonstrated from (60) to (63), with a simple distributive property, it is possible to obtain a direct relationship between the reduced and the full Schur complement. These direct relations cannot be made in a model with two or more fields to be reduced. The choice of one reduced field, associated with the Schur complement, permits a MOR based on Galerkin projection in a similar way of processing of an optimization problem without constraints. The benefits of such an approach are the reduction of matrix operations and the application of several existing adaptive strategy already developed. Back to the formulation, the reduced system can be developed as (38) and (39), ie, where Δz Ψ x represents the approximate value of increment of the special variable "z" from the range space of the primal variable. In the same way, an approximate value of increment of the dual variable from the range space of the primal variable is obtained as (40), ie, At this point, it is clear that These variables are approximations that depend on the quality of the representation from the Galerkin projection of the primal variable space x . A more detailed explanation related to the error associated with both types of variables is explained in the next section. The next subsection is related to the dual approach of the Galerkin projection.

MOR of the dual field
In this section, the main hypothesis (a reduction of one field) considers the decomposition of the Lagrangian multiplier (dual variable), ie, where is the local range space operator related to the dual space, and is the local null space operator. The main variables are arranged in the following decomposition: When these global Galerkin operator and the global null space operator are replaced in (52), the parcel of the error is represented by the term As in the previous subsection, the term of (72) is extremely useful to measure the error related to the full model, including the application in an adaptive strategy to reduces the error of the MOR. Finally, the terms related to the projection are Applying the Schur complement as developing in (42) to (44), where Again, with a simple distributive property from (76) and (77), it is possible to obtain a direct relationship between the reduced and full Schur complement. As in the previous subsection, these direct relations cannot be made in a model with two or more fields to be reduced. Once again, the choice of one reduced field leads to a similar way of processing of an optimization problem without constraints. In the specific case of this work, the explicit convex approximation functions generate diagonal Hessians; this matrix format reduces, even more, the number of operations to reach the inverse matrices. Returning to the formulation development, the reduced system can be developed as (45) and (46), ie, where Δz represents the approximate value of increment of the special variable "z" from the range space of the dual variable. In the same way, an approximate value of increment of the primal variable from the range space of the dual variable is obtained as (47), ie, As in the previous subsection, at this point, it is clear that These entire variables are approximations that depend on the quality of the representation from the Galerkin projection of the dual variable space . A more detailed explanation related to the error associated with both types of variables is described in the next section.

The IPM algorithm with MOR and the KKT condition
Until now, as described in the first subsection of this work, the MOR occurs in the linearization of the manifold, the vector space  u . Such approach is common to computational mechanics; optimization without constraint (using the classical Newton-Raphson) or optimization with equality constraints can apply snapshots from previous linearized spaces without major concerns. However, in optimization methods with inequality constraints, the KKT condition is a factor that difficult the use of previous snapshots to develop a Galerkin projection, with low error, to the next step.
The source of such problem is the barrier parameter; the level of variation (that is an empirical choice) can change abruptly the values of the gradient and Hessian. In the next section, more precisely in the subsection related to the sensibility analysis, it is demonstrated that even small variations of the barrier parameter can result in meaningful variations to the primal or dual variables. This work adopts the increment of the primal and dual variables as snapshots; because of that, to avoid the increase of the error in the projection, the option adopted is to restart the Galerkin projection for each new value of barrier parameter. In Flowchart 1, the IPM algorithm with MOR is described.
The additional variables in the flowchart are The difference between the full and MOR algorithms resides in obtaining the new increment through the development of the Newton Method. The application of the MOR in Flowchart 1 justifies the snapshots acquired from the increments. It is important to highlight that the option of discard the snapshots with a new barrier parameter is based in some possibilities of high variations demonstrated in the sensibility analysis from Section 4; it is not possible to ensure that this option always result in a lower number of corrections to the Galerkin projection.

Error related to the primal variable
In the previous section, the parcel of the error related to the projection from the primal variable is obtained from the vector described in (56); based on this vector, it is possible to associate the error of the Galerkin projection related to the primal variable. To achieve that, the idea is to isolate the vector from the null space of the primal variable. This approach is possible, exposing all the equations from (52) with null space of the primal variable different of a null vector, ie, Isolating the projection value in (84), Now, replacing the projection value (86) in (85) and isolating by field, Related to the error analysis, the operators A 12(x) , A 13(x) , A 14(x) , and a 1(x) are less relevant; the description of such operators can be seen in Appendix B. The relevant operator is related to the vector of the null space, and in this case, it is the operator A 11(x) , where In the null space perspective, it is possible to see that all the other terms must belong to In addition, if the Galerkin projection converges to the full space, This occurs because of the parcel of (88), ie, Based on this parcel, it is possible to realize that the error is dependent of the Hessian of the Lagrangian function related to the primal variable, and the curvature of this function in the primal field must be convex to ensure the convergence of the error (the matrix must be invertible). This Hessian, described in (22), is partially dependent on the dual variables, and consequently has some levels of dependency to fulfill the KKT condition.

Error related to the dual variable
Furthermore, in the previous section, the parcel of the error related to the projection from the dual variable is obtained from the vector described in (72); based on this vector, it is possible to associate the error of the Galerkin projection related to the dual variable. As in the previous subsection, the idea is to isolate the vector from the null space of the dual variable. This approach is possible, exposing all the equations from (52) with null space of the dual variable different of a null vector, ie, Isolating the projection value in (94), Now, replacing the projection value (96) in (92) and isolating by field, As in the previous subsection, all the operators that are not related of the null space vector from the dual variable are described in Appendix B. The relevant operator is related to the vector of the null space, and in the case of (97), it is the operator A 13( ) , where In the same way, related to the others operators, it is possible to achieve replacing the projection value (96) in (93) and isolating again the fields, ie, and the operator A 23( ) , associated with the null space, Finally, replacing the projection value (96) in (95) and isolating the fields, Results in the operator A 33( ) , ie, As in the previous subsection, it is possible to see that all the other terms must belong to the null space operators In all operators, when the Galerkin projection converges to the full space, This occurs because of the parcel that is common to all operators of the null space, ie, Moreover, in this case, it possible to take two conclusions from the error related of Galerkin projection in the dual variable. The first conclusion is that the null space is directly dependent of the variables that compound the KKT condition. The direct dependency of the dual and slack variables results, in the case of IPM algorithm, a higher dependency of the barrier parameter value. The consequence of such dependency is demonstrated in the last subsection. In the next subsection, it is possible to notice the advantage of the Galerkin projection applied in one field per time when compared with the application of both Galerkin projections at the same time.

Error related to the primal and dual variable at the same time
In the last two subsections, the analysis of the null space operators demonstrates that a Galerkin projection in a unique field generates no direct error between the primal and dual variables. To understand better the advantage of the selection between the primal and dual fields, once per time, it is important to develop an error analysis with the Galerkin projection in two fields at the same time. To start such analysis, the fields that are reduced are defined, ie, Now, the global Galerkin projection and the global null space operator are defined in (106), they are replaced in (52); the parcel of the error is represented by the term Like in the previous two subsections, the idea is to isolate the two vectors from the null space of the primal and dual variables. Again, this approach is possible, exposing all the equations from (52) with null spaces of the primal and dual variables different of null vectors, ie, Developing a similar manipulation described in the previous two subsections, but now isolating the two null space and the other two fields in one equation, The questions related in how to achieve Equation (112), as such as the values of all operators, can be answered in Appendix B. Once again, the important operators to the error analysis are related to both null spaces It is possible to see that all the others terms must to belong to When the range space converges to the full space, It is possible to realize that (113) and (114) contain the parcels described in (91) and (105), respectively. In the point of view of the dual variable, the null space operator, even with both projections applied at the same time, remains with the same behaviors described in the previous subsection. However, in the point of view of the primal variable, the behavior of the operator is different when compared to the first subsection. Now, the operator in (113) accumulates both errors related to the Galerkin projection in the primal and dual variables. In the case that an increase of error in the dual field exists, both fields are affected. Furthermore, in the case that an increase of error in the primal field exists, such field must deal with the error from the dual field too, resulting more time processing to reduce the error of the range space in both cases.
As can be seen in the first subsections, when the Galerkin projection is applied in the primal field, the error in the primal field is only partially dependent on the KKT condition. These property permits, when the barrier parameter is not properly chosen, lower errors in the primal field when compared to the dual field, which this last one is completely dependent on the barrier parameter. The level of dependency with the KKT condition of each parcel, described in (91) and (105), can be seen in more details in the next subsection, through the sensibility analysis.

IPM algorithm sensibility that affects the MOR error
In this subsection, the idea is a sensibility analysis related to the operators that compound the parcels (91) and (105). Based on such parcels, this subsection analyzes the sensibility from the Hessian of the Lagrangian function related to the primal field, the Lagrange multipliers (dual field) and the slack variables. Although the range space converging to the full space eliminates all null space operators, before this happens, the operation between the inverse of the reduced operator with the original operator can amplify, in this meantime, the distance between the range and full spaces.
Others operators from the parcels (91) and (105) are the Galerkin projections from the primal and dual variables. Related to these operators, there are two points to consider. The first point is that the sensibility analysis depends on the adaptive strategy from the MOR; the choice by the generation or update is made by eigenvectors, singular vectors, or Ritz vectors, as such as the criteria of selection change the sensibility analysis. The second point is that the sensibility of primal and dual variables changes the snapshots that generate the Galerkin projections, becoming the process extremely complicate to analyze. Because of these both particularities of the adaptive strategy, this subsection focuses on the operators from the infeasible primal-dual IPM with artificial variables.
The sensibility analysis starts in the IPM algorithm. In such algorithm, as described in Flowchart 1, the barrier parameter must achieve the value lower than a default tolerance; in the meantime, the barrier parameter assumes different values, decreasing until the default tolerance value. Considering now two sequential barrier parameters, from step "j" to step "j + 1", as the imposition of two incomplete KKT conditions, Assuming the difference between (118) and (117) results It is possible to realize that the sensibility (variation of the dual variable and slack variable) of the IPM algorithm is dependent on the variation (or decrease) from the barrier parameter. Usually, there are adaptive strategies to determine the variation of the barrier parameter to achieve a better convergence; however, there is no exact function to describe the barrier parameter variation. Because of that, the convergence of the IPM algorithm, as such as the error described by the parcel represented by (105), is dependent on an empirical estimation of the variation from the barrier parameter. In addition, it is important to highlight that a low difference between barrier parameters is not able to ensure a low variation of both variables; when low differences are imposed, three situations can occur, ie, Based in (120) and (121), it is possible to realize that the low variation can occur in only one variable. Although the term from (105) is compound with the inverse of the dual variable, the sensibility analysis represented by (120), (121), and (122) are enough to demonstrate the dependency of (105) related to the barrier parameter. In the case of the term from (91), the dependency of the barrier parameter is not as clear as (105), but it is possible to realize that some terms that compound the Hessian of the Lagrangian function related to the primal variable have direct influence with the barrier parameter.
To make the sensibility analysis clearer, the terms related to the superior and inferior limit of the primal variable from (22) are represented by Replacing (123) in (22), As the previous analysis of sensibility, considering now two sequential barrier parameters, from step "j" to step "j + 1", and replacing (117) and (118) for each step, respectively, Again, assuming the difference between (126) and (125) results where Splitting Equation (127) in one parcel referent to the variation from the primal variable, and another one related to the variation from the dual and slack variables, where Operator (130) is directly dependent on the variation from the primal variable, and it is clear that, when the variation of the primal variable tends to zero, such operator vanishes. However, in operator (131), the analysis of the variation follows the conditions from (120), (121), and (122), ie, In the same way, a low difference between barrier parameters cannot ensure a higher convergence related to the primal variable. Another condition, with no relation to the barrier parameter, is the high variation generated by (123). Such condition comes from the definition of the problem, which is in two extremes situations that can be achieved from a specific value of the primal variable, ie, The variations described in (135) and (136) are related to the constraints of the problem. In this case, it is possible to obtain a low variation from (131) only if condition (134) is fulfilled. Based on what has been shown so far, it is possible to realize that a selection of the barrier parameter variation that ensures a fast convergence is essential to achieve quickly the low variation of parcels (91) and (105). Once a low variation is achieved, the snapshots from the previous points (vector spaces) from the manifold  are closer to each other, ensuring Galerkin projections with a better representation of the current vector space  u .
To reduce the MOR error and achieve a faster convergence, this work recommends the application of two empirical approaches. The first one is, obviously, an adaptive strategy to the barrier parameter variation that leads to a convergence closer in the quadratic rate. The second approaches is the normalization of the objective and constraints functions; in case of both types of functions have similar proportion, the variation of the Lagrange multiplier (dual variable) is lower, resulting also in a lower variation of the term (131). In consequence of different behaviors from both kind of null space operator and by the lack of objective approach to control them, the next section of this work suggests an adaptive strategy that changes the field of Galerkin projection, according to the comparison between errors.

Introduction to the MOR and adaptive strategy
Although the MOR presented until now can be applied as a surrogate model for several steps, the proposal mentioned previously is to generate an initial guess to the PCG, in a similar way as described in the works of Tromeur-Dervout and Vassilevski 10 and Markovinovic and Jansen. 15 However, different of the works cited previously, there is a selection based on error related to each projection. Although such an approach can be applied in several types of analysis, an efficient processing time only is achieved to a low-cost inverse Hessian. An option to obtain a low-cost inverse Hessian is to apply an explicit convex approximation function as described in the work of Bruyneel et al. 29 Such approximation function generates a diagonal (explicit property) and invertible (convex property) Hessian.
The purpose of this section is to use the models developed in Sections 2 and 3 and, through an accurate initial guess, to provide an improvement of convergence to the PCG method. To facilitate the understanding of how to select the models, this section is divided into three subsections. The first subsection is how to generate and restart (if necessary) the operator of the Galerkin projection. The second subsection is related on how to select the initial guess from a Galerkin projection with lower error. The last subsection is how to update and correct the operator of the Galerkin projection. Flowchart 3 and Flowchart 4 show, respectively, the MOR with adaptive strategy and only the MOR inside the previous flowchart.

Initialization and restart of the Galerkin projection
To develop the initial projection operators, it is necessary to collect snapshots of both fields of the problem. In this work, where all the analyses are developed online, the gathering starts after the acquisition of the new value of the barrier parameter. In a specific barrier parameter value, a predetermined number of iterations k v is selected; each iteration obtains two snapshots, ie, the first one from the increment of the primal and the second one from dual variable. The result of the procedure described in Algorithm 2 is the initial Galerkin projection for each field (commands inside box means that the process can be parallelized). After such procedure, the MOR is initialized. During the MOR, the matrices described in Algorithm 2 suffer updates; these updates store the last k v snapshots. Such matrices are applied when a Galerkin projection operator reaches a maximum dimension w max n or w max m ; then, a new classical POD is applied to obtain a new Galerkin projection. In Algorithm 3, this procedure is demonstrated more clearly. Between these two processes described in Algorithm 2 and Algorithm 3, the Galerkin projection operators suffer variations of the dimensions to reach the next points (vector spaces  u ) in the manifold . To obtain the variation of dimension, the solution by the projection of the fields (range space) and the full solution of the fields are necessary. The next subsection describes how to obtain both solutions.

Selection of the initial guess and PCG
In this work, the solution of the reduced space provides the initial guess of the PCG. A contribution of this research is the initial guess selected by the Galerkin projection with lower error. Before the description of the process, it is important to highlight that the efficiency of this approach depends on inverse Hessian with a simple solution, as described previously. Consider that, in such approach, both reduced Schur complements are processed per iteration, as described in Algorithm 4.
Once the range of space vectors of the primal and dual variables is obtained, the original dimension of each field is analyzed. The original dimension of the field selects the Schur complement applied in the process. As a recommended criterion, the lower original dimension selects if the Schur complement is applied to the primal or dual variable. Based on the previous criterion, if the Schur complement with a primal variable is chosen (n < m), the initial guess is between Equations (65) and (81). In (65), the equation comes from the range space of the primal variable, and Equation (81) comes from the range space of the dual variable.
However, if the Schur complement with a dual variable is chosen (n > m), the initial guess is between Equations (67) and (79). In (67), the equation comes from the range space of the dual variable, and Equation (79) comes from the range space of the primal variable. This work tests two kinds of criteria; the criteria to select between Equations (65) and (81), or between Equations (67) and (79). Based on that, to the problem solved with the Schur complement from the primal variable, To the Schur complement from the dual variable, In the first criterion, the values from (137) and (138) are related to the lower error from the initial guess or the lower error from the Galerkin projection, ie, Moreover, the second criterion, which comes from (56) and (72), is defined by the minimum Galerkin error Unlike (139), the selection from (140) does not depend on the field related to the Schur complement to solve the full system. In Algorithm 5, the procedure of the selection described previously is demonstrated more clearly.
In this stage, with a solution in the range space and the full solution, it is possible to find the solution of the null space; such space is applied to update and correct the Galerkin projection operator to both fields. The next subsection describes such a process.

Update and correction of the Galerkin projection
As described in the previous subsection, with the solution from the range and full space, it is possible to obtain the solution of the null space by the difference of both. To update the Galerkin projection, it is applied the concept of Ritz vector. In this case, the Ritz vector is developed by the orthonormalization (by Gram-Schmidt process) of the null space vector, ie, In relations (141), the orthonormalization occurs related to the Galerkin projection operator vectors. Consequently, after the development of the Ritz vectors, such vectors are added to the Galerkin projection. In the Galerkin projection operators, the Ritz vector increases the number of columns; the number of columns represents the dimension of the reduced space. During the update of the Galerkin projection, the norm of errors is calculated, and the error can be developed at the same time of the update from the Galerkin projection operators. This procedure is demonstrated more clearly in Algorithm 6.
The process represented by Algorithm 6 is repeated per iteration. During the IPM algorithm, the repetition of this process, described in Flowchart 3, is possible to generate a nonnull kernel in the reduced matrices from Equations (64) and (78), resulting in singular matrices. To avoid this behavior, the determinants of the matrices from Equations (64) and (78) are calculated; when a determinant reaches a value extremely close to zero, the last Ritz vector is discarded. In this case, the Galerkin projection operator has the number of columns reduced in one dimension with the attempt to restore the bijectivity of the matrices cited previously. The reduction of the number of columns, and consequently the dimension of the reduced space, is repeated until the values of the determinant from each matrix increases again. This procedure is demonstrated more clearly in Algorithm 7.
All the processes described in this section, from Flowchart 2 to Algorithm 7, show that this work is concerned to simplify the adaptive strategy. The classical POD and Ritz vectors are basic strategies to the Galerkin projection generation, and the reduction of this one is developed only to ensure the solvability of the reduced systems. These works assume this approach to focus the analysis in the reduction of one field per time and the improvement of the convergence from the PCG by the selection of the initial guess generated by the MOR. The next section shows the results of simulations based on this adaptive strategy.

Definition of the examples
As mentioned in the last section, this work developed two examples to demonstrate the advantages and disadvantages of the adaptive strategy described previously. These examples are an extension of the example described in the work of Svanberg 26 ; the main difference is related to the number of constraint functions, but properties as convexity and limited domain in the primal variable are the same, ie, A modification is made related to the aforementioned work, 26 the explicit convex function adopted is the classical method of moving asymptotes (MMA), as described in the work of Svanberg 25 and in the work of Bruyneel et al. 29 It is important to highlight also, different of the works of Wächter and Biegler 24 and Herskovits, 30 the IPM algorithm is processing without line-search. These both modifications are made to facilitate the reproducibility of the simulations. More details can be seen in Appendix C and Appendix D.
Moreover, it is important to mention that all reduced linear systems are solved by Cholesky decomposition, and the preconditioned of the CG method, different of more advanced approaches as in the works of Ghosh et al, 31 Risler and Rey, 32 and Anzt et al 33 is always the diagonal of the Schur complement matrix. Related to the simulated models, this work compares four kinds of simulations: the simulation with the full model, MOR with a fixed field, MOR with field selection based of initial guess, and MOR with field selection based of the Galerkin error. Each simulation measured the time and the number of iterations of the PCG per step.
Related to the procedure, both examples are simulated following the procedures.
1. Both examples are simulated by applying the Schur complement in the primal form (59) and dual form (75). The application of the Schur complement occurs in separate simulations and they are indicated by the term "primal" or "dual" in the title of the graph. 2. In the full model, when simulated with the Schur complement in the primal form, it applies the previous value of Δx as initial guess; when simulated with the Schur complement in the dual form, it applies the previous value of Δ as initial guess. 3. The simple MOR simulation is made with the Galerkin projection always in the same field. When the MOR is simulated with the Schur complement in the primal form, it applies the value of (65) as initial guess; when simulated with the Schur complement in the dual form, it applies the value of (79) as initial guess. The result from this simple MOR simulation is described in the legend of the graphs by "x" to the primal form and "lambda" to the dual form; 4. The MOR simulation proposed in this work is made with Galerkin projection alternating the fields. When the MOR is simulated with the Schur complement in the primal form, it can apply between the equations from (65) and (81) as initial guesses that are dependents of a criterion; when simulated with the Schur complement in the dual form, it can apply between the equations from (67) and (79) as initial guesses that are dependents of the same criterion. The result from this MOR simulation is described in the legend of the graphs by "x" and "lambda". 5. In the MOR simulation proposed in this work, the criteria to select the initial guess can be made between (139) and (140). The first one is based on the minimum error between the initial guess and the final result of the PCG; the second is based on the minimum error between the Galerkin error fields. In the graphs, result from (139) is described in the legend by "initial guess error", and the result from (140) is described in the legend by "Galerkin error". More details are presented in Algorithm 5. 6. Related to the demonstration of the results in the graphs, one part of graphs is demonstrated related to the steps; a step represents the output of Flowchart 1, and in the graph, it is represented by 7. Another part of them is demonstrated related to the internal iterations of the IPM algorithm that occur inside the algorithm from Flowchart 1, and in the graph, it is represented by Value of the algorithm ] k,i n total = number of steps that exist the iteration k.
The internal iterations result in a step; each step has a specific number of internal iterations. Because of that, the number of iterations showed in the graph is determined by the step with more iteration. In the graphs, the differentiation of both approaches is showed through the name of the horizontal axis as "step" or "iteration". The different representations help to understand global aspects of the results (steps), as well as local aspects, which occur inside the IPM algorithm (iterations).

Example 1
The first example is a simulation of a problem with more design variables (primal) than constraints (dual). Such example represents the majority of problems in computational mechanic problems with multifields. The data are In (145), they are respectively the number of design variables, the number of constraints, the initial value of each design variable, and the maximum dimension of each reduced space.

Example 2
The second example is a simulation of a problem with more constraints (dual) than design variables (primal). Such example does not represent the majority of problems in computational mechanic problems with multifields, but it has a behavior that is important to show for this work. The data are In (145), they are respectively the number of design variables, the number of constraints, the initial value of each design variable, and the maximum dimension of each reduced space.

Results of the examples
The first results in Graph 1 show the total number of iterations in the PCG per step. It is possible to realize, in these examples, that Graphs 1B and 1D show reductions related to the MOR with the selection of the fields if compared to the full model and MOR only with the dual field. However, in Graphs 1A and 1C, the reductions are only meaningful related to the full model; the selection of the field does not show relevant reduction during the simulation. The explanation of the different behavior by the different application of the Schur complement is described at the end of these examples through an analysis of error. Table 1 shows the accumulative number of iterations in the PCG.
For now, a more direct analysis is made related to the processing time; the second type of results compare the processing time of the PCG from the full model with the processing time of the PCG from the three others reduced models. However, for a fair comparison, each reduced model is added the processing time of the reduced linear system and the norm of errors.
Graphs 2A, 2B, 2C, and 2D do not show relevant changing of results when compared to the PCG iteration in Graphs 1A, 1B, 1C, and 1D. This means, at least for the dimensions selected to the both reduced space, that the additional processing time caused to the adaptive strategy is not relevant to the total processing time of the problem. Table 2 shows the accumulative processing time during the entire simulation.   Table 3 shows the accumulative processing time, only of the adaptive strategy, during the entire simulation. Although the selection of the fields increases the processing time, it is realized that the adaptive strategy described in this work shows some better results compared to the full model and the MOR with only one reduced space; because of that, the next results are focused on the MOR with the selection of the field. The next type of results shows the number that each field is selected per step based on the adaptive strategy described in the last section.
In Graphs 3A and 3C, it is possible to notice that the selected field of the norm from initial guess error changes considerably (the curves intersect) if compared to the selected field of the norm from the Galerkin error. This behavior occurs in the primal field because the adaptive strategy imposes a reduction described in Algorithm 7. The reduction in the local      Galerkin projection acts directly of the norm from initial guess error. Such behavior occurs only in some parcels of the norm from the Galerkin error; the result is a norm much less sensitive than the norm from initial guess error. Table 4 shows the total value of the selected fields from the examples.
To explain better what happens to the primal field, the next type of results demonstrates the reduction suffered by the adaptive strategy described in Algorithm 7; they show the average of the reduced dimension of the MOR in each step.
The decrement of reduced dimensions, related on to initial guess error, in Graphs 4A and 4C, changes completely the selection of the fields described in Graphs 3A and 3C. However, the adaptation to the new local Galerkin projections does not change considerably the number of the PCG iterations, as demonstrated in Graphs 1A and 1C. Table 5 shows the total average of the reduced dimension from the examples.
In the next type of results, the norm from initial guess error is not analyzed; the choice is made by the major instability, gives by the correction in Algorithm 6, and low changing from the previous results. The next type of results shows the evolution of the selected fields during the internal iterations of the IPM algorithm; the criterion of the selection is the Galerkin error, as described in (140).
Analyzing Graphs 5A, 5B, 5C, and 5D, it is possible to realize that the error related to the primal field is lower during the most part of the simulation. In the middle of the simulation, it is the only part that the error of the dual field becomes Reduced dimension -average Ex.   lower; such behavior can be explained, first, by the variation of the Galerkin operator and, second, by the variation of the barrier parameter during the IPM algorithm. Related to the barrier parameter, the next type of results shows the evolution during the internal iterations of the IPM algorithm. It is possible to realize, in Graphs 6A, 6B, 6C, and 6D, that the value of 1E-3 results in a faster convergence (the number of this value is reduced related to the value of 1 and 1E-6) and consequently reduces the error of the dual field. However, the value of 1E-6 in the final part of the simulation is critical in this case; in the IPM algorithm, a small value of the barrier parameter, at the end of the simulation, increases the condition number of the system. Moreover, an ill-conditioned system induces a poor convergence to the PCG, and a poor initial guess (higher Galerkin error-damages, even more, the situation. Table 6 shows the accumulated selection of the field according to the barrier parameter. Finally, the last type of results shows the PCG iterations per internal iteration of the IPM algorithm: As expected, as the final part of the simulation, as shown in Graphs 7A, 7B, 7C, and 7D, the values of the PCG iterations increase considerably. Furthermore, it is possible to notice that the number of the PCG iterations in the dual examples reduces dramatically when the primal field is selected, specially, in the last half of simulation (Graph 5B against Graph 7B  and Graph 5D against Graph 7D). Related to the primal examples, no relevant reduction of PCG iterations is noticed during the simulations, not even in the middle of the simulation (the lowest error of the dual field); an explanation is that the barrier parameter is not small enough to result in an ill-conditioned system, and consequently, a better initial guess is not so relevant to reduce the number of PCG iterations. Until now, at least in these two examples, the different behaviors from both Schur complements are clear. The part of the explanation can be found in Section 4, through the error and sensibility analysis. Subsection 4.1, specifically in (91), shows that the error of the primal field is not directly dependent on the KKT condition, suffering a less impact of the barrier parameter selection, unlike the dual field, showed in subsection 4.2, specifically in (105), which a direct dependency of the error with KKT condition is verified. The influence of the KKT condition is even more evident at the end of simulation, when low values of barrier parameters cause a slower convergence; because of that, no relevant reduction in the PCG iteration is noticed in the first half of the simulation.
Nonetheless, it is important to highlight in these two examples that, despite an error analysis and a sensibility analysis, they are not enough to ensure that the selection of the dual field always results in higher error and lower convergence; however, it is possible to conclude that there exist higher risk of high error and low convergence in the dual field, if compared to the primal field, because of the direct dependency of the empirical selection from the barrier parameter.

CONCLUSION
In this work, it possible to realize that the combination of the Galerkin projection in one field, applied to the infeasible primal-dual interior point Method, can be simplified with the approach of the Schur complement. However, both Galerkin projections still have, even in different levels, an error dependency of the KKT conditions. Nevertheless, a MOR using Galerkin projection can be applied with a considerable improvement of the processing time if some approaches are made.
1. As demonstrated in both examples, in problems with explicit approximation convex functions, the adaptive strategy using the selection of the field can generate a considerable reduction in the processing time in a simulation with the Schur complement in the dual form. This approach has a broad application in topology optimization, especially in problems with a great number of inequality constraints. In the case of the topology optimization, the use of explicit approximation convex functions result in diagonal Hessians from the objective function; the calculus of the inverse hessian to apply in the Shur complement in the dual form becomes extremely cheap. 2. In problems with only differentiable convex objective functions, where the diagonal Hessian is not possible, an interesting approach can be the Galerkin projection only in the primal field. Although the simulation of this kind of problem is not showed in this work, all the formulations presented here, including the error and sensibility analysis, permit such approach. In the Schur complement only in the primal field, it does not require the calculation of the inverse Hessians from the objective function; the only inverse that is required comes from the KKT condition that is naturally diagonal. Moreover, as demonstrated in Section 4, the Galerkin error from the primal field is less dependent from the barrier parameter selection. Because of these arguments, the processing time from the MOR based in Galerkin projection only in the primal field can result in some reduction of time, if compared to the full model. Such an approach can be applied in some kind of problems in mechanics, as contact mechanics and plastic behavior of the material.
It is important to mention that the analysis of the MOR described in this work is much wider than what is simulated here. An analysis with different adaptive strategies to the barrier parameter variations, to improve the convergence of the IPM algorithm and reduces the MOR error, could generate another work. Furthermore, problems with minimum maximum (min-max), reduction of other fields (related to "y", for example), and the combination of reduced fields, which the parcel of error is related to the primal and dual field as demonstrated, are only some of the possibilities to develop in the next research works. Finally, to further research works, problems with original convex function, as described in the last topic, are interesting applications to the IPM algorithm with Galerkin projection.

APPENDIX A
The Schur complement of a matrix divided into four blocks can be obtained in two ways: to isolate x 2 , the multiplication of the first row by A 21 A 11 −1 and subtracting the same related to the second row are necessary; and to isolate x 1 , the multiplication of the second row by A 12 A 22 −1 and subtracting the same related to the first row are necessary. These manipulations result in ∀A 11 ∈ R n 1 x n 1 , A 22 ∈ R n 2 x n 2 , A 12 , A 21 T ∈ R n 1 x n 2 , b 1 , x 1 ∈ R n 1 , b 2 , x 2 ∈ R n 2 . (A1) To simplify the number of variables, it is defined Now, considering the equivalent system as (30)  ] [ x 1 ∀ a 13 , a 31 T ∈ R n 1 x 1 , a 23 , a 32 T ∈ R n 2 x 1 , a 33 , x 3 , b 3 ∈ R. (A3) The system described in (31) can be achieved in two parts. In the first part, the first row and the second row from (A3) being separate and manipulated, ie, Arranging b 1 to isolate x 3 and replacing the terms by the terms from (36) and (34), Replacing also A 1 , x 1 , and x 3 by the equivalent term from (32) and replaces in (A4) results in the first row of (31), ie, Related to the second row, the second row and the third row from (A3) being separated and manipulated, ie, Arranging b 1 to isolate x 1 and replacing the terms by the terms from (37) and (35), Replacing also A 1 , x 1 , and x 3 by the equivalent term from (33) and replaces in (A7) results in the second row of (31), ie, Now, considering the equivalent system as (41) defined in Section 2.4, A 23 , A 32 T ∈ R n 2 x n 3 , A 33 ∈ R n 3 x n 3 , a 34 ∈ R n 3 x 1 , x 3 , b 3 ∈ R n 3 , x 4 ∈ R. (A10) Now, multiplying the first row by A 31 A 11 −1 , the second row by A 32 A 22 −1 , and subtract both related to the third row, Replacing also A 3 , a 34 , x 3 , and x 4 by the equivalent term from (43) and (44), and replaces in (A11), results in the first row of (42), ie,

APPENDIX B
The rest of the operators from Section 4.1 are In Section 4.2, the rest of the operators, from the first row, are In the second row, Moreover, in the third row, Related to Section 4.3, now is shown the development to achieve the two null space operators. The development starts isolating the projection value of the primal variable in (108), ie, In addition, isolating the projection value of the dual variable in (109), Replacing (B6) in (B5), it leads to Now, isolating the projection value of the dual variable in (110), ie, Replacing (B8) in (B7), where the barrier parameter variation is +1 = 0.001 .
Furthermore, the normalization of the functions is defined as In the next step, the problems start with the same values (only x starts with the value of the previous step). The values of the matrices W, P, and Q are These matrix values are the same as the classical example defined in the work of Svanberg. 26 The values from the matrices R k and T k defined in this work are