Rapid Multi-Objective Optimization of Periodically Operated Processes Based on the Computer-Aided Nonlinear Frequency Response Method

The dynamic optimization of promising forced periodic processes has always been limited by time-consuming and expensive numerical calculations. The Nonlinear Frequency Response (NFR) method removes these limitations by providing excellent estimates of any process performance criteria of interest. Recently, the NFR method evolved to the computer-aided NFR method (cNFR) through a user-friendly software application for the automatic derivation of the functions necessary to estimate process improvement. By combining the cNFR method with standard multi-objective optimization (MOO) techniques, we developed a unique cNFR–MOO methodology for the optimization of periodic operations in the frequency domain. Since the objective functions are defined with entirely algebraic expressions, the dynamic optimization of forced periodic operations is extraordinarily fast. All optimization parameters, i.e., the steady-state point and the forcing parameters (frequency, amplitudes, and phase difference), are determined rapidly in one step. This gives the ability to find an optimal periodic operation around a sub-optimal steady-state point. The cNFR–MOO methodology was applied to two examples and is shown as an efficient and powerful tool for finding the best forced periodic operation. In both examples, the cNFR–MOO methodology gave conditions that could greatly enhance a process that is normally operated in a steady state.


Dynamic Process Intensification and Forced Periodic Operation Analysis
The idea of intensifying a process by switching from a steady-state (SS) to a forced periodic (FP) operation was first introduced by Douglas and Rippin in 1966 [1], and it was further expanded by Douglas [2], Horn and Lin in 1967 [3], and Bailey and Horn in 1971 [4]. Though in the beginning, the investigations of FP operations were limited to theoretical studies, in the next 25 years, there were at least 29 different experimental dynamic process intensification (PI) investigations of industrially important processes [5]. The initial investigations studied the FP operations with the modulation of one input variable: the inlet flow rate or the chemical composition (reactant concentration). The studies focused on enhancing the reaction rate, selectivity, and product distribution [5]. However, except for the reverse flow reactors, no remarkable improvements, i.e., process enhancements, were achieved by applying PI concepts in the time-domain. This was the case until the emergence of sophisticated process 1.
A systematic framework for the identification of possible dynamically intensified processes is required.

2.
Various mathematical and numerical challenges while optimizing and computing the cyclic processes need to be addressed. 3.
New techniques that encompass the previous points should be part of intuitive software tools that can be used both in academia and industry.
The idea of this research was to address the beforementioned challenges using a completely new approach based on an analytical PSE tool for predicting the time-average behavior of periodic processes-the Nonlinear Frequency Response (NFR) method [35]. Recently, this method was upgraded by developing a new PSE software application, leading to the so-called computer-aided Nonlinear Frequency Response (cNFR) method [36].

The NFR Method for Evaluating the Potential of Process Intensification through Forced Periodic Operations
The Nonlinear Frequency Response (NFR) method is an analytical PSE tool based on the concept of higher-order frequency response functions (FRFs) [37]. The NFR method allows users to transfer a dynamic mathematical model of a weakly nonlinear system from the time domain to the frequency domain. In the frequency domain, the system is mathematically represented with a set of the first, second, and higher-order FRFs [37], which are functions of frequency. Though an infinite sequence of 3 of 20 FRFs is needed for full representation of the system, it was shown that good approximations can be achieved by using only the first and second-order FRFs [38].
For an easier understanding of the NFR method, Figure 1 shows a schematic representation of a chemical reactor with forced, sinusoidal modulation of an input x around its SS value x s . For a stable system, after going through a transient, any output from the reactor will eventually reach a quasi-(periodic)-steady-state y qs . Most processes in chemical engineering, including chemical reactors, are weakly nonlinear (all nonlinear terms in their model equations are continuous and differentiable, i.e., can be represented in the Taylor series form). The NFR method can only be applied to stable, weakly nonlinear systems that can be represented with convergent Volterra series. Another limitation is that the analyzed systems must not have multiple steady states. For such systems, the quasi-steadystate of the response is obtained as a sum of its SS value, the first harmonic (of the same frequency as the input), a theoretically indefinite number of higher harmonics, and a nonperiodic term (the so-called DC component), as shown in Figure 1. The DC component represents the difference between the time-average of the output in the FP regime and its SS value and is the measure of PI owing to the FP operation.
stable system, after going through a transient, any output from the reactor will eventually reach a quasi-(periodic)-steady-state . Most processes in chemical engineering, including chemical reactors, are weakly nonlinear (all nonlinear terms in their model equations are continuous and differentiable, i.e., can be represented in the Taylor series form). The NFR method can only be applied to stable, weakly nonlinear systems that can be represented with convergent Volterra series. Another limitation is that the analyzed systems must not have multiple steady states. For such systems, the quasi-steady-state of the response is obtained as a sum of its SS value, the first harmonic (of the same frequency as the input), a theoretically indefinite number of higher harmonics, and a nonperiodic term (the so-called DC component), as shown in Figure 1. The DC component represents the difference between the time-average of the output in the FP regime and its SS value and is the measure of PI owing to the FP operation.
The NFR method based on the concept of higher-order frequency response functions (FRFs) [39] enables the fast and easy estimation of the DC component. For the sinusoidal modulation of a single input , with amplitude and frequency , the DC component can be approximately evaluated by using only the asymmetrical second-order FRF , , relating the output and input : The sign of , , which is the same as the sign of , , , answers whether the FP operation would result in process improvement or not, while its value gives a very good estimate of the extent of the improvement. When two input signals ( and ) are modulated, with amplitudes and and a phase angle between them, then: the DC component can be approximately calculated as a sum of contributions of the inputs and , separately, and the contribution of their cross-effect [40]: The NFR method based on the concept of higher-order frequency response functions (FRFs) [39] enables the fast and easy estimation of the DC component. For the sinusoidal modulation of a single input x, with amplitude A and frequency ω, the DC component can be approximately evaluated by using only the asymmetrical second-order FRF G (2) x,y (ω, −ω) relating the output y and input x: The sign of DC x,y , which is the same as the sign of G x,y (ω, −ω), answers whether the FP operation would result in process improvement or not, while its value gives a very good estimate of the extent of the improvement.
When two input signals (x 1 and x 2 ) are modulated, with amplitudes A 1 and A 2 and a phase angle ϕ between them, then: the DC component can be approximately calculated as a sum of contributions of the inputs x 1 and x 2 , separately, and the contribution of their cross-effect [40]: x2,y (ω, −ω)+ 2· x1x2,y (ω, −ω) + sin ϕ·Im G (2) x1x2,y (ω, −ω) where G (2) x1x2,y (ω, −ω) is the cross asymmetrical second-order FRF. From Equation (4), it follows that the simultaneous modulation of two inputs has a higher PI potential in FP operations than single-input modulations (Equation (1)), as the sign of the cross-term can be easily adjusted by correctly selecting ϕ [40]. Due to this, by modulating the system with two inputs with an optimal phase difference, improvements become possible even in cases when modulating one or both inputs, separately, leads to process deterioration [41].
It is important to note that different performance criteria of interest (e.g., conversion, yield, selectivity, and productivity) of FP operation can be directly calculated from the time-average values of the outputs obtained using the NFR method. Using Equations (1) and (4), they can be expressed in the algebraic form as functions of the input modulating parameters (frequency, amplitudes, and phase shift).
One of the most important steps of applying the NFR method is the derivation of the needed FRFs. Though the procedure for the analytical derivation of FRFs is well-established and has been used and presented in many publications [35,[39][40][41][43][44][45][46][47][48][49][50][51][52][53][54], it is rather complex, time-consuming, and requires certain mathematical skills in the complex domain. As such, the NFR method can be unappealing to beginners, especially when applied to more complex systems.
The recently developed computer-aided Nonlinear Frequency Response (cNFR) method [36] is an upgrade of the NFR method using a software application for the automatic derivation of the FRFs of interest. It has an easy-to-use interface where the user can define the dynamic model equations of the system and automatically generate MATLAB files for the desired FRFs and DC components [36]. The generated MATLAB files contain codes with algebraic expressions of the derived FRFs and DC components in the vector form, and they enable very fast simulation of the key system performance criteria and rapid rigorous global and local optimization. No time-consuming and numerically expensive methods are required for solving the systems of differential equations. In this way, the cNFR method provides answers to three key questions about the potential of the investigated FP operation:

1.
Which input or combination of inputs should be periodically modulated? 2.
What are the optimal forcing parameters of the modulated inputs: frequency, amplitudes, or phase shift (ω, A 1 and A 2 , or ϕ, respectively)? 3.
What is the expected process improvement if we switch from SS to FP operation (a maximal increase in the performance criteria of interest)?
By giving answers to these questions in a fast and easy manner, the cNFR method overcomes most challenges listed by Baldea and Edgar. As is shown in this paper, the cNFR method has the potential to become a very strong software PSE tool for analyzing the dynamic PI potential.
The cNRF software tool was explained in detail by Živković et al. [36]. Thus far, it was used for the analysis of the FP operation of an isothermal continuous stirred tank reactor (CSTR) with the simultaneous modulation of feed concentration and flow-rate [24], for experimental identification [55], and for the evaluation of the FP operation of an electrochemical process [41].
In this work, a new methodology for the rigorous multi-objective optimization (MOO) of forced periodic (FP) operations is presented. It is based on the cNFR method and cost-benefit indicator analysis.

The cNFR-Multi-Objective Optimization Methodology for Periodically Operated Processes
In Figure 2, the cNFR-MOO methodology is shown schematically. It consists of steps done by the user (gray rectangles) and stages that are a direct result of the application of the cNFR software tool (white rectangles). Explanations of how the dynamic model equations, modulated inputs, and outputs of interest are defined is given in a publication by Živković et al. [36], where the cNFR tool was presented. After using the cNFR tool, in a matter of minutes, the user gets automatically generated MATLAB files that contain all FRFs, DC components of interest, and a system characteristic equation that can be used for the determination of system stability ( Figure 2). the user (gray rectangles) and stages that are a direct result of the application of the cNFR software tool (white rectangles). Explanations of how the dynamic model equations, modulated inputs, and outputs of interest are defined is given in a publication by Živković et al. [36], where the cNFR tool was presented. After using the cNFR tool, in a matter of minutes, the user gets automatically generated MATLAB files that contain all FRFs, DC components of interest, and a system characteristic equation that can be used for the determination of system stability ( Figure 2).
System stability is an important prerequisite for the use of the cNFR method [36] and should be used as one of the algebraic constraints in MOO. The system is stable as long as all roots of the characteristic equation derived by the cNFR method remain negative [36]. Therefore, this necessary stability constraint ( Figure 2) is already provided as a function to the user and should be included in the MOO code. Other constraints can also be defined if needed. The objective functions reflect important system criteria (e.g., reactor capacity, productivity, and selectivity) and can differ from case to case. In general, a cost indicator (CI) is set as a combination of the DC components of system variables that are undesired and expensive from the process point of view. On the contrary, the benefit indicator (BI) is defined as a combination of the DC components of the most beneficial variables in the system. An objective function can also contain a combination of both BIs and CIs. MOO constraints and objective functions expressed through CIs and BIs are defined in the MATLAB code by the user (Figure 2 [25]. It is important to point out that in the cNFR-MOO methodology, we optimize the steady-state (SS) and forcing parameters all in one step. This is different than in classical approaches of optimizing System stability is an important prerequisite for the use of the cNFR method [36] and should be used as one of the algebraic constraints in MOO. The system is stable as long as all roots of the characteristic equation derived by the cNFR method remain negative [36]. Therefore, this necessary stability constraint ( Figure 2) is already provided as a function to the user and should be included in the MOO code. Other constraints can also be defined if needed.
The objective functions reflect important system criteria (e.g., reactor capacity, productivity, and selectivity) and can differ from case to case. In general, a cost indicator (CI) is set as a combination of the DC components of system variables that are undesired and expensive from the process point of view. On the contrary, the benefit indicator (BI) is defined as a combination of the DC components of the most beneficial variables in the system. An objective function can also contain a combination of both BIs and CIs. MOO constraints and objective functions expressed through CIs and BIs are defined in the MATLAB code by the user (Figure 2). The main reason that MOO is preferred over SOO is that it enables users to examine multiple optimal solutions for different values of competing objective functions of interest. Additionally, it provides results in cases where SOO can lead to mathematically infinite ( lim (CI1·CI2), etc.) [25].
It is important to point out that in the cNFR-MOO methodology, we optimize the steady-state (SS) and forcing parameters all in one step. This is different than in classical approaches of optimizing periodic operations, which are usually performed in two steps: first by finding the optimal SS point and then by finding the optimal forcing parameters (frequency, amplitudes, and phase shift) of the input modulations around that SS. The reason for simultaneous optimization of SS and forcing parameters is that the DC component-defining the process improvement-depends on the SS values of the process variables, and it is possible that perturbing the system around a sub-optimal instead of an optimal steady state can result in higher improvements. Thus, it is possible to find a periodic operation that gives the best possible results regarding the chosen performance criteria. As a result of such optimization, a complete set of parameters is obtained: the SS point and the forcing parameters of the modulated inputs.
The final step shown in Figure 2 consists of results analysis and additional FRF and DC component simulations if the further evaluation and understanding of the FP operation are needed. The cNFR method is applied to MOO in Section 3 on two examples: (1) an isothermal continuous stirred tank reactor (CSTR) with a simple reaction mechanism and (2) an electrochemical reaction (ECR) process.

Problem Formulation
In Example 1, the following irreversible n-th-order reaction is taking place in an isothermal continually stirred tank reactor (CSTR) with a constant volume, where A is the reactant, P is the product, and k r is the rate constant of a simple n-th order reaction rate: For an isothermal CSTR, only material balance equations are needed to fully mathematically describe the system: where t is time, V is the volume of the reaction mixture, F is the volumetric flow rate, C Ai is the reactant inlet concentration, and C A and C P are the outlet concentrations of reactant and product, respectively. The molar flow rates for reactant inlet and product outlet can be expressed as: This simple case was already used as an example for the analysis of potential PI by FP operations, first by using the standard NFR method [39,40] and then by the cNFR method [36]. Here, this example is used to illustrate MOO based on the cNFR method for finding the optimal periodic operation for the case of simultaneous sinusoidal modulations of the feed concentration C Ai and flow rate F: where index s denotes the variable steady-state (SS) and A 1 , A 2 , and ϕ are forcing parameters. For the purposes of optimization, dimensional frequency, ω, is used instead of the nondimensional, ω nd , according to the following equation: The optimal SS operation is used as a base case for comparison.

Formulation of the Objective Functions and Constraints
The two objective functions that are used are the maximized product yield Y P and the reactor volumetric capacity κ. For periodic operations with variable flow rates, these values are defined based on the time-average values of the inlet and outlet molar flow rates [43]: By expressing the time-average outlet molar flow rate of the product, based on the corresponding DC component DC FC Ai ,MF P , Equation (14) becomes [36]: where the subscript s denotes the SS values of the variables. Product yield (Equations (14) and (16)) is the first benefit indicator (BI1). The second benefit and also a cost indicator (BI2/CI) is the reactor volumetric capacity (Equation (15)), and its value is maximized to get a smaller reactor volume, as a cost indication, and higher capacity, as a benefit indication.
In the current analysis, the only constraint is that the system must be stable. As explained in our previous publications [36,40], the condition of the reactor stability can be expressed in the form of the following inequality: − C A,s ·F s − C A,s n ·V·k r ·n < 0 (17)

Optimization Variables and Criteria Formulation
In Table 1, the optimization variables for both the SS and FP operations are listed. The SS CSTR operation has only SS optimization variables, i.e., volumetric flow rate (F s ), inlet concentration (C Ai,s ), and the reaction mixture volume (V). On the other hand, the FP operation has additional dynamic optimization variables: the amplitudes of the input modulations (A 1 and A 2 ), the frequency of modulation (ω), and the phase angle between modulated inputs (ϕ). The bounds are set so that all variables are optimized between approximately 0 and 1, and then they are scaled by multiplying with the corresponding scaling factor, as shown in Table 1. Therefore, C Ai,s is optimized between approximately 0 and 20 kmol/m 3 , V between 1 and 100 m 3 , ω between 0 and 2π rad/s (corresponding to 1 Hz), and ϕ between 0 and 2π rad. Higher forcing frequencies are not considered because it is not realistic for flow rate or concentration change to be faster than 1 Hz.
The MOOs were performed using MATLAB 2019b (MathWorks, Natick, MA, USA) and a nondominated sorting genetic algorithm (NSGA-II). The criteria for the NSGA-II were set according to the values listed in Table 2. To obtain the optimal SS values, the system of nonlinear algebraic equations was solved with the trust-region dogleg algorithm, with the function termination tolerance (TolFun) set to 1 × 10 −12 . The total number of optimization variables differed between the SS and FP operations ( Table 2) and corresponded to the variables listed in Table 1.
After executing both MOOs on a device with a 1.9 GHz Intel i7-8665U 4-core processor, an Intel 620 UHD graphics card, and 32 GB RAM, the total time needed for the SS operation optimization was 142 s and 189 s for the FP operation. Therefore, the dynamic optimization took less than a minute longer than the SS optimization even though the population was set to 500 and the tolerance was set to 1 × 10 −4 ( Table 2).

Results and Discussion
The optimization results are shown in Figure 3 in the form of Pareto fronts. It is desirable to have both Y P and κ as high as possible.
The MOOs were performed using MATLAB 2019b (MathWorks, Natick, MA, USA) and a nondominated sorting genetic algorithm (NSGA-II). The criteria for the NSGA-II were set according to the values listed in Table 2. To obtain the optimal SS values, the system of nonlinear algebraic equations was solved with the trust-region dogleg algorithm, with the function termination tolerance (TolFun) set to 1 × 10 −12 . The total number of optimization variables differed between the SS and FP operations ( Table 2) and corresponded to the variables listed in Table 1. After executing both MOOs on a device with a 1.9 GHz Intel i7-8665U 4-core processor, an Intel 620 UHD graphics card, and 32 GB RAM, the total time needed for the SS operation optimization was 142 s and 189 s for the FP operation. Therefore, the dynamic optimization took less than a minute longer than the SS optimization even though the population was set to 500 and the tolerance was set to 1 × 10 −4 ( Table 2).

Results and Discussion
The optimization results are shown in Figure 3 in the form of Pareto fronts. It is desirable to have both and as high as possible.
(a) (b) Two different process enhancement zones, marked with roman numerals in Figure 4a, could be identified by analyzing the optimal solutions (colored stars and gray triangles): • Zone I (blue stars) was in the region of very low κ and high Y P , i.e., κ ≈ 0 and Y P > 90%. In Figure 4a, Zone I is shown with a dashed rectangle. Both SS and FP operations showed the same trend of achieving the highest product yields Y P , at the lowest CSTR volumetric capacities, κ (Figure 4a). This was expected as the highest Y P could be achieved with negligible flow rates and in the smallest reactors, i.e., with the highest residence time (Figure 4a). Zone I showed that as Y P increases, the influence of PI through input modulation, or system nonlinearity enhancement, potential decreased when compared to the benefits of a higher reactor residence time. The result was complete reactant conversion for the SS and FP operations but with minimal production potential (low κ). Zone I was of no interest in this research because it lied in the region of negligible κ and did not enhance the process significantly if the FP operation was utilized.

•
In Zone II (red stars, Y P < 90%, Figure 4a), the FP operation showed PI potential while outperforming the optimal SS operation. Since Zone II was of interest for the dynamic enhancement of the CSTR process, it is also shown, in part, in Figure 4b. The same Pareto front was displayed on a different scale in Figure 4b with chosen cases for analysis: the SS case and the three FP cases of FP1, FP2, and FP3 (denoted with small circles in Figure 4b). For the chosen cases, optimal Pareto results can be seen in Table 3. As visible from Figure 3 and Table 3, the FP operation gave significantly higher Y P and κ when compared to the optimal SS operation. The optimal operating parameters were rapidly determined thanks to the cNFR-supplied DC function files. In Appendix A, a detailed analysis of the MOO can be found. In the next example, the cNFR-MOO methodology is used to find an optimal FP operation for enhancing electrochemical oxygen reduction reaction, which is of industrial relevance.

Problem Formulation
The oxygen reduction reaction (ORR) is utilized in fuel cells, batteries, and some electrolysis processes. It is also of strong industrial relevance in processes that involve chlorine production [55]. The kinetics of the ORR were previously studied using the cNFR method by Kandaswamy et al. [55]. The first-and second-order symmetrical FRFs were verified for 0.1 M NaOH solution [55]. Živković et al. explained in detail how the FRFs were derived using the cNFR method [36]. In a subsequent study, the cNFR method was used to evaluate electrochemical process improvement if FP operation is utilized [56].
The ORR is taking place on the surface of a silver rotating disc electrode, with the assumed single-step mechanism [55]: In total, six equations describe a single-step mechanism shown in Equation (18): 1.
The electrical charge balance: where c dl is the double-layer capacitance, η the electrode overpotential, F the Faraday's constant, r the ORR reaction rate expression, and curr the current density, which is the main output of interest.

2.
The potential balance: where E is the electrode potential (which is a possible modulated input), E θ is the standard electrode potential, and R Ω is the ohmic resistance.

3.
The reaction rate of the ORR: where k app is the apparent rate constant of the ORR according to the equation: In Equation (22), k e1 is the ORR kinetic constant, α is the transfer coefficient, R is the universal gas constant, and T is the reaction temperature, which is assumed to be constant.

4.
The intermediary oxygen concentration derived from the discretized mass balance Equation [55]:

5.
The electrode boundary layer thickness: where ω r is the electrode rotation rate, D is the oxygen diffusivity in the boundary layer, and ν is the kinematic viscosity of the NaOH solution.
In the additional investigation, the second-order asymmetrical FRF and the corresponding DC component of the current density, curr (when electrode potential, E, is periodically modulated), were experimentally verified by Živković et al. [56]. In the same research, the verified DC component was used to show that there can be significant ORR process improvement when E is periodically modulated input, x, around its SS value, E s (according to equation in Figure 1) [56]: Here, the evaluation of the ORR FP operation was extended further by conducting rigorous cNFR-MOO.

Objective Functions, and Constraints Formulation
In a preliminary evaluation, the FP ORR operation was considered to be superior to the SS operation if the absolute mean value of the current density, curr , was higher and the absolute mean value of the overpotential, η , was lower [56]. Therefore, CI as η , and BI as curr were seen to be two objective functions. These objective functions could be evaluated based on the corresponding asymmetrical second-order FRFs using the following equations [56]: In Equations (26) and (27), G E,η (ω, −ω) and G (2) E,curr (ω, −ω) are the asymmetrical second-order FRFs for the periodically modulated potential E (input variable), and η and curr are outputs of interest, respectively. The cNFR software gives the user MATLAB function files with the derived FRFs and the respective DC components in an algebraic, vector form. Thus, η and curr can be computed in an easy and fast manner. It should be noted that FRF G (2) E,curr (ω, −ω) was experimentally verified in a study by Živković et al. when higher current densities were achieved by using forced periodic operation [56].
One constraint for both the SS and FP operations is that the system must be stable. Therefore, the roots of the characteristic equation obtained by the cNFR equation have to be negative. The characteristic equation was derived in the same manner as for the previous example, and, due to its numerous and long equation terms, is not shown here. The MOO of the FP operation also has two additional constraints that satisfy the following inequality: Equation (28) was introduced to ensure that at any time the potential at the electrode surface could not be higher than 1 V or lower than 0.1 V because, at those conditions, no ORR can physically occur.

Optimization Variables and Criteria Formulation
The optimization variables for both SS and FP operations are given in Table 4. Two electrode operating variables were optimized in both operations: the electrode surface SS potential, E s , and electrode rotation rate, ω r . E s cannot be lower than 0.1 or higher than 1 V, thus satisfying the inequality interval constraint (Equation (28)). The minimum possible ω r was set to 400 rpm, and the maximum was set to 2500 rpm, well-within the technical limitations of the rotating device described by Kandaswamy et al. [55]. ω r was scaled to be optimized between 0.16 and 1, and its real value was calculated by multiplying with the maximal allowed value of 2500 rpm, as shown in brackets in Table 4.
The FP operation has two additional, dynamic optimization variables: the amplitude of FP modulation A, and the frequency of modulation, ω. A was optimized between 0 and 1, while frequency was optimized between 10 −3 and 10 3 Hz, with its exponent as the optimization variable. The variable was further multiplied by 2π to obtain the radial frequency, ω, in rad/s (Table 4).
The MOOs of the SS and FP operations were performed in MATLAB 2019b by using the same algorithm (NSGA-II). The criteria for NSGA-II, which are listed in Table 5, were the same as in Example 1 ( Table 2), except that the number of optimized variables was now different (two for SS operation and four for FP operation) and parallel computing was used because of more constraint evaluations.
After conducting MOO on the same device as in Example 1, the total MOO time for the SS operation was 118 s and 324 s for the FP operation. The dynamic optimization took less than four minutes longer than the SS optimization.

Results and Discussion
The results of the MOOs for the SS and FP operations are shown in Figure 4. The SS operation (gray triangles in Figure 4a) gave an increased BI at a higher CI, as expected. The FP operation (colored stars in Figure 4a) showed a similar trend with a sudden discontinuity for the CI between approximately 0.5 and 0.65 V, i.e., at a BI of approximately 61 A/m 2 . operation and four for FP operation) and parallel computing was used because of more constraint evaluations. After conducting MOO on the same device as in Example 1, the total MOO time for the SS operation was 118 s and 324 s for the FP operation. The dynamic optimization took less than four minutes longer than the SS optimization.

Results and Discussion
The results of the MOOs for the SS and FP operations are shown in Figure 4. The SS operation (gray triangles in Figure 4a) gave an increased BI at a higher CI, as expected. The FP operation (colored stars in Figure 4a) showed a similar trend with a sudden discontinuity for the CI between approximately 0.5 and 0.65 V, i.e., at a BI of approximately 61 A/m 2 .   Three distinctive PI zones, marked with different colored stars and surrounded with dashed rectangles and Roman numerals in Figure 4a, could be identified by analyzing the Pareto front: • Zone I (blue stars) was in the region of very low CI and BI, i.e., η < 0.35 V and curr < 4 A/m 2 . In Zone I, both operations gave similar BI and CI and produced almost no current density. This zone was marked by the dominant kinetic regime and negligible ORR in the system.

•
In Zone II (red stars), for 0.35 < η < 0.50 V and 4 < curr < 61 A/m 2 , or until the discontinuity, the FP operation showed a strong PI potential and outperformed the SS operation as CI increased. This zone was of the highest interest in the dynamic enhancement of the ORR process and is also shown, in part, in Figure 4b.

•
Zone III (green stars), which was in the range of high CI and BI, i.e., η > 0.65 V and curr > 61 A/m 2 , showed no ORR enhancement as both SS and FP operations gave the same BI and CI. In this zone, incremental improvements in BI were very expensive as they could only be achieved at exponential increases in CI.
In Table 6, optimal Pareto results for selected cases marked with circles in Figure 4b, along with the optimized operational values, are given. As can be seen both in Figure 4 and Table 6, the optimal FP operation outperformed the optimal SS operation in Zone II, both in BI and CI values. It should be noted that if MOO was not used but SOO with a combined CI and BI (in form BI/CI) was, the optimal solution would be in Zone III, thus resulting in no improvement with FP operation. By using MOO, the user can see all possible PI enhancement zones and identify operating FP conditions for which a 334% higher current density was achieved when compared to the SS case. A more detailed analysis of the Figure 4 and Table 6 results can be found in Appendix B. The discontinuity line that occurred at curr ≈ 61 A/m 2 , and 0.5 V < η < 0.65 V, was void of any optimal FP solutions ( Figure 4a). This can be explained by the sudden change of sign of the respective FRF functions (Equations (26) and (27)). A more detailed analysis of the discontinuity and Zones I, II, and III can be found in Appendix C.

Conclusions
In this paper, we introduced a completely new approach to optimizing the forced periodic operations of processes that are classically operated in a steady-state regime. This was also the first time that the simultaneous optimization of both forced periodic and its corresponding operational steady-state point was performed for a chemical process in the frequency domain. The main contribution of this approach was based on combining multi-objective optimization techniques with the Nonlinear Frequency Response (NFR) method, which enables direct estimation of time-average values of the outputs of periodically operated processes. Thanks to the NFR method, the objective functions could be defined as algebraic expressions. These expressions depend on the steady-state point around which the system is perturbed, on one hand, and the forcing parameters of the input perturbations on the other. Recently, the NFR method was upgraded to the computer-aided Nonlinear Frequency Response (cNFR) method using a software tool for the automatic derivation of the frequency response functions (FRFs), which are crucial for applying the NFR method. The main advantages of the cNRF-based multi-objective optimization of periodic processes can be summarized as the following: The objective functions are defined in the form of algebraic expressions defining the time-average behavior of the periodic process, so there is no need for numerical solutions of the nonlinear dynamic model equations. Consequently, the computing times for the optimization of the periodic and the steady-state operations are of the same order of magnitude.

2.
Due to the automatic derivation of the FRFs using the cNRF software, defining and evaluating the objective functions of interest are fast and easy tasks. Periodic operations with one or two modulated inputs can be essentially treated in the same way.

3.
The new approach performs the optimization of the steady-state point and the forcing parameters in one step. In this way, in some cases, it is possible to find a periodic operation around a sub-optimal steady state that would be superior, not only to the steady-state operation but also to any periodic operation around the previously established optimal steady state.
The methodology for cNFR-assisted cost-benefit indicator multi-objective optimization was clearly defined and illustrated in two different cases: an isothermal continual stirred tank reactor (CSTR) with a simple reaction mechanism and two simultaneously modulated inputs-the reactant concentration and flow rate of the feed stream-and a process with electrochemical oxygen reduction reaction (ECR), with the modulation of only one input-the electrode potential. As shown in both CSTR and ECR examples, forced periodic operation can greatly intensify the process under the right operating conditions. When two inputs were modulated, the determination of the optimal phase angle, which was possible due to the cNFR-generated function files, was of the utmost importance. The product yield, the primary CSTR benefit indicator, was doubled under the FP operation when compared to the optimal SS operation with the same volumetric reactor capacity. Similarly, in the previously experimentally verified ECR example, the PI of the optimized FP operation led to a current density 3.34 times higher than in the optimized SS operation. The usefulness of the cNFR as a new PSE tool did not end there, as it enabled the detailed analysis of the obtained Pareto fronts through easy and quick simulations of the key FRFs. With the cNFR method, it was possible to identify all PI zones and then focus on optimal solutions with the highest potential for process improvement. In both examples, the speed of the multi-objective optimizations was remarkable (the dynamic optimizations were finished in several minutes), even at tightened tolerances and increased genetic algorithm populations.
As shown, the methodology for cNFR-based multi-objective optimization, presented in this work, can be a powerful PSE tool that could bring a new impulse to the promotion and implementation of forced periodic operations in the chemical industry. The FP operation with the simultaneous modulation of F and C Ai enabled a CSTR with a doubled volumetric capacity for the same product yield when compared to SS operation.
When the FP1 case (red circle in Figure 4b) is compared to the SS case (gray circle in Figure 4b) for the same product yield (approximately 61%), the increase of volumetric reactor capacity in the FP operation can be seen to have been approximately 100%.

2.
Likewise, the FP operation have higher product yield when operated in the reactor with the same volumetric capacity as the SS operation.
If we compare the FP3 case (red circle in Figure 4b) and the SS case (gray circle in Figure 4b) for the same volumetric capacity (approximately 0.1 kmol/m 3 /min), the product yield for FP operation can be seen to have been approximately 18% higher. This difference became even bigger for lower product yields.

3.
In FP operation, it is possible to find operational values that give both a higher product yield (more benefit) and an increased reactor capacity (fewer costs) than an optimal SS operation.
The SS case (gray circle in Figure 4b) have a lesser product yield and allowed for a smaller reactor capacity when compared to the FP2 case (red circle in Figure 4b). Thus, the FP operation can be superior to SS operation when both BI and CI are compared.
From the data presented in Table 3, the following additional findings can be noted:

4.
Even if SS and FP operations are carried out with similar residence times, the FP operation can still be superior to the SS operation in BI and CI.
If we compare the FP2 and SS cases (both shown in Figure 4b, with values listed in Table 3), Y P and κ can be seen to have been higher in the FP operation, while the residence time was almost the same (162 min for the SS case and 163 min for the FP2 case). However, under the FP operation, the residence time fluctuated with the modulated volumetric flow rate from the lowest value of 84 min to the highest value of 45 h. The key in using the enhancement potential of the fluctuating residence time is in making sure that the phase angle between modulated inputs is carefully selected (optimized value).

5.
The optimal phase difference between the two modulated inputs in all FP cases was around 6 rad.
The modulation of F was slightly behind the modulation of C Ai (by approximately 16 • ) for all selected FP operation cases listed in Table 3. The optimally selected phase difference between modulated inlets led to superior results when the FP was compared to SS operation. A correctly selected phase angle enhanced the reaction in the CSTR and was the reason why FP1 had a better CI, FP3 had a better BI, and FP2 had better both BI and CI when compared to the SS case (Table 3).

1.
The FP operation could give the same current density as the SS operation at a reduced CI (higher overpotential), and this reduction enlarged as more current density was produced.
If the FP1 case is compared to the SS case (Figure 4b), it can be seen that the reduction in costs was approximately 18%. This reduction became bigger as solutions for higher BIs were compared, up until Zone III ( η > 0.65 V).

2.
The FP operation allowed for a dramatically increased BI at the same CI values.
The FP3 case produced approximately 230% more current density when compared to the SS case ( Figure 4b) while at the same overpotential. As can be seen in Figure 4a, this trend was increased at higher CIs and held until the discontinuity that happened at approximately 0.5 V.

3.
Similar to Example 1, it was possible to find optimal parameters of the FP operation that have both higher BI and lower CI values when compared to the optimal SS operation.
When the FP2 case is compared to the SS case in Figure 4b, it can be seen that the FP operation produced approximately 110% more current density at 11% less CI (higher overpotential).
Additional findings are listed after analyzing the optimal Pareto results for selected cases given in Table 6 and marked with circles in Figure 4b For both the SS and FP operations, the optimal electrode rotation rate was at its maximum allowed value.
As already explained by Živković et al. [56], a higher ω r leads to a smaller diffusion layer thickness at the electrode surface and a faster ORR. This, in return, gives a higher BI.

5.
More intensified FP operation cases have lower SS values of electrode potential and higher amplitudes of input modulation.
The most intensified FP operation case, FP3, had the lowest E s , followed by FP2 and then FP1. The finding that smaller values of E s could lead to an increased BI seems counterintuitive. Živković et al. found that the highest ORR process improvement can be expected in a kinetically controlled region, i.e., an increased E s leads to an increased curr [56]. However, that was concluded for a significantly lower ω r (1600 rpm), so it can be assumed that the mass transport was slowing down the ORR for the higher E s more than in the optimized cases shown in Table 6. Higher amplitudes for more intensified cases were expected because the upper bound for E was the same for all cases (1 V). Lower E s values, in return, allowed for more expressed perturbations with higher A values (Equation (28)).

6.
The frequency of modulation for intensified cases was slightly increased but of the same order.
Cases FP1, FP2, and FP3 had a frequency of the same order, although it was optimized between 10 −3 and 10 3 Hz ( Table 4). The value of frequency was slightly increased for more intensified cases but well-within the bounds predicted by our previous FP operation evaluation study, where it was found that PI can happen only for frequencies lower than 1 kHz, with the highest PI potential between approximately 10 and 100 Hz [56].

Appendix C. Analysis of PI Zones and Discontinuity in Example 2: ECR
Additional investigation was done because of discontinuity between Zone II and III described in Section 3.2.4 (shown in Figure 4a). Mathematically, the change of shape in Pareto front (Figure 4a), which was visible in transition from Zone I to II and then again from Zone II to III, was the consequence of the sign change of the DC component used to calculate the absolute time-average current density (BI), i.e., DC E,curr (Equation (27)). As already explained in Section 1.2, the shape and sign change of a DC component will be the same as its corresponding FRF if only one input is modulated (Equation (1)). In this example, the FRF of interest was the second-order asymmetrical FRF for the current density, E,curr (ω, −ω), as shown in Equation (27). If the sign of G (2) E,curr (ω, −ω) is positive, there will be PI potential and FP will outperform the SS operation. If the FRF value is zero, there will be no PI potential, and FP and SS operations will give the same process efficiency. If FRF has a negative sign, the SS operation will be superior and no PI by this input FP modulation will be possible. In Figure A1, the ORR enhancement potential was studied for three different modulation frequencies by simulating the respective FRF, , , , for different values of the SS potential, (between 0.2 and 0.9 V). The results are shown in Figure A1 as a function of absolute time-average overpotential, i.e., CI.
As can be seen in Figure A1a, the red (10 Hz) and blue lines (10 mHz modulation) have a positive value up until CI ≈ 0.56 V, when their value becomes zero. This means that for a CI between approximately 0.3 and 0.56 V (Figure A1a), the system will outperform itself in the SS operation if it Figure A1. The oxygen reduction reaction (ORR) enhancement potential as a function of time-average overpotential (CI) at three different modulation frequencies (10 mHz as blue, 10 Hz as red, and 10 kHz as a green line): (a) the second-order asymmetrical frequency response function (FRF) for current density and (b) the second-order asymmetrical FRF for overpotential.
In Figure A1, the ORR enhancement potential was studied for three different modulation frequencies by simulating the respective FRF, G  Figure A1 as a function of absolute time-average overpotential, i.e., CI.
As can be seen in Figure A1a, the red (10 Hz) and blue lines (10 mHz modulation) have a positive value up until CI ≈ 0.56 V, when their value becomes zero. This means that for a CI between approximately 0.3 and 0.56 V (Figure A1a), the system will outperform itself in the SS operation if it is forced to periodically operate (FP operation) around the same SS. This should not be confused with the comparison of FP and SS operations in Figure 4a, which had different optimal operating SS values. From Figure A1a, it is clear that the performance of the FP operation will be better at 10 Hz than at 10 mHz because the red line will give slightly higher G (2) E,curr (ω, −ω) values. For a CI between approximately 0.56 and 0.8 V, (Figure A1a), the sign of G (2) E,curr (ω, −ω) for lower frequencies (the red and blue lines) will be negative, while for high frequencies (the green line), the sign will be zero. This means than no PI will be possible and if FP operation is used; it should be done so at higher frequencies, but the result will still be the same as in the corresponding SS. At CI ≈ 0.62, the red and blue lines for G (2) E,curr (ω, −ω) will reach a minimum. The numbers extracted and analyzed from Figure A1a correspond to findings of the discontinuity line in Figure 4a, which was placed in the CI range between 0.50 and 0.65 V.
For higher CI values, G E,curr (ω, −ω) will be zero and no PI can be expected if FP operation is used. The same thing can be concluded if the second-order asymmetrical FRF for overpotential, G E,η (ω, −ω) is simulated ( Figure A1b). The sign change of G (2) E,η (ω, −ω) also happens at a CI of approximately 0.56 V, same as with G (2) E,curr (ω, −ω) in Figure A1a. In Figure A2, all Pareto solutions for the optimal frequency of modulation are plotted against the CI ( Figure A2a) and BI ( Figure A2b). Depending on the zone they belong to (Figure 4a), the solutions are presented in Figure A2 with blue (Zone I), red (Zone II), and green (Zone III) stars. Between Zone II and III in Figure A2a, there is a clear gap that corresponds to the discontinuity line shown in Pareto front (Figure 4a) at 61 A/m 2 (boundary line between Zone II and III in Figure A2b).
All zones show different optimal frequency ranges. Zone I ( Figure A2), which was characterized by a negligible ORR and a dominant mass transfer, had an optimal frequency below 10 Hz. The most interesting zone for analyzing PI by FP operation, Zone II, had an optimal frequency in the range between 10 and 50 Hz ( Figure A2). This followed the previous evaluation study on the PI of the ORR [56]. As can be seen in Figure A2, Zone III, which offered no PI under the FP operation (Figure 4a), had a significantly higher optimal frequency than Zone I and II, up to 1 kHz (which was the upper bound for frequency; see Table 4). This followed Figure A1, where it was shown that for a higher CI, the only way to use the FP operation is at higher frequencies but with no process improvement In Figure A2, all Pareto solutions for the optimal frequency of modulation are plotted against the CI ( Figure A2a) and BI ( Figure A2b). Depending on the zone they belong to (Figure 4a), the solutions are presented in Figure A2 with blue (Zone I), red (Zone II), and green (Zone III) stars. Between Zone II and III in Figure A2a, there is a clear gap that corresponds to the discontinuity line shown in Pareto front (Figure 4a) at 61 A/m 2 (boundary line between Zone II and III in Figure A2b).
All zones show different optimal frequency ranges. Zone I ( Figure A2), which was characterized by a negligible ORR and a dominant mass transfer, had an optimal frequency below 10 Hz. The most interesting zone for analyzing PI by FP operation, Zone II, had an optimal frequency in the range between 10 and 50 Hz ( Figure A2). This followed the previous evaluation study on the PI of the ORR [56]. As can be seen in Figure A2, Zone III, which offered no PI under the FP operation (Figure 4a), had a significantly higher optimal frequency than Zone I and II, up to 1 kHz (which was the upper bound for frequency; see Table 4). This followed Figure A1, where it was shown that for a higher CI, the only way to use the FP operation is at higher frequencies but with no process improvement achieved when compared to SS operation.