# Introduction Singular Optimal Control Problem (SOCP) consists in determining the control variable profiles that minimize an objective function, subject to algebraic and differential constraints. In the last decade, a significant increase of control techniques in the industrial context was observed. The reason for this is mainly due to the high popularity of dynamic simulation tools and the existence of a competitive global market, in which environmental constraints and demanding market specifications require a continuous improvement of process operation. Dynamic optimization enables an automatic decision-making procedure. Therefore, as it gets established as an useful and trustworthy technology, other industrial applications are driven forward even more efficiently, such as: the addressing of hard constrained problems, the synthesis of chemical reactors networks, the uncertainties description in multiple period problems and the development of tools such as automatic differentiation (Biegler et al., 2002). In order to solve this kind of problem, several numerical methods have been proposed (Bryson and Ho, 1975). They are usually classified according to three broad categories, regarding their underlying formulation: direct optimization methods, Pontryagin's Minimum Principle (PMP) based methods, and HJB-based methods. The PMP approach is based on the optimal control theory and requires the numerical solution of multipoint boundary value problems involving state and adjoint (costate) variables. The main difficulty associated with using this type of method is the initial estimate for the costate variables (Costa, 1996;Biegler et al., 2002). In the context of chemical engineering, a typical example of a SOCP is the fermentation process, where the substrate concentration can be maintained at a fairly low level and unfavorable effects of a high concentration, such as growth inhibition, can be avoided. This phenomenon leads to unimodal reaction rate expressions, which exhibit a maximum point with respect to a single reactant concentration or in terms of two or more reactant concentrations. Although only one single control variable, in the form of the feed rate, may appear to characterize a simple optimal control problem, considerable difficulties have been reported in the determination of the optimal feed rate policy for fedbatch processes, due to the intrinsic nonlinearity of these systems (Hong, 1986;Modak et al., 1986;Modak and Lim, 1989;Fu and Barford, 1993;Xiong and Zhang, 2003). In this problem, the usual objective considered in the optimization of a fed-batch bioreactor is to maximize the metabolite production or the yield, that is, the production per unit of substrate fed (Hong, 1986). Traditionally, during engineering system design, the model, the vector of design variables, and the parameter vector are considered free of errors, i.e., they do not contain uncertainties. However, more realistically, small variations in the vector of design variables may cause significant modifications in the vector of objective functions (Ritto et al., 2008). As a consequence, the system to be optimized can be very sensitive to small changes in the vector of design variables, and thus, small variations in this vector can cause significant changes in the vector of objective functions (Ritto et al., 2008). In this context, it is important to determine a methodology that produces solutions less sensitive to small variations in the vector of design variables. Solutions with this characteristic are called robust solutions and the procedure to find these solutions is named Robust Optimization (Taguchi, 1984). In this contribution, the Multi-objective Optimization Differential Evolution (MODE) algorithm (Lobato, 2008), associated with the Effective Mean Concept-EMC (Deb and Gupta, 2006) is applied to determine the feed substrate concentration in fed-batch penicillin fermentation process. This robust multiobjective singular optimal control problem consists of maximizing the productivity and minimizing the operation total time. In the post-processing stage of the results, the criterion adopted to choose a point on the Pareto curve is the overall profit. This work is organized as follows. Sections 2 and 3 presents the mathematical description of the SOCP and the mathematical model that describes the fed-batch penicillin fermentation process, respectively. Section 4 shows a brief review about the MODE algorithm. The EMC strategy considered to deal with uncertainties is presented in Section 5. The results obtained are presented in Section 6. Finally, the conclusions are outlined in Section 7. # II. # Optimal Control Problem The solution of an OCP consists in the determination of the control variables profiles that maximize or minimize a measure of performance. The OCP performance index is given by: (1) where ? and L are the first and second terms of the performance index, respectively. The objective is subject to the implicit Differential-Algebraic Equations (DAE) system: (2) with initial conditions assumed consistent and given by: ( A comparison among methods for solving the OCP had great attention around the first part of the eighties with the development of numerical methods, appropriate to a more restricted class of problems, identified mainly by the differential index (Brenan et al., 1996). The indirect strategy for solving the OCP is based on variational principles. These conditions, from the Pontryagin's Minimum Principle (Bryson and Ho, 1975), generate a set of Euler-Lagrange equations, which are boundary value problems (BVPs), inherently formed by the DAE, regardless of whether the problem is restricted or not. Some difficulties in the OCP solution must be highlighted: (i) the existence of end-point conditions or region constraints implies in multipliers and associated complementary conditions that significantly increase the difficulty of solving the BVP by the indirect method; (ii) the existence of constraints in the state variables and the application of the slack variables method may produce DAE of higher indexes, regardless of the constraint activation status, even in problems where the number of inequality constraints is equal to the number of control variables; and (iii) the Lagrange multipliers may be very sensitive to the initial conditions. The direct approach, on the other hand, uses the control parameterization (sequential method) or the state and control parameterizations (simultaneous method), transforming the original problem into a finite dimensional optimization problem. By all means, the implementation of direct methods is simpler because it does not demand the generation of the costate equations, which, at very least, duplicates the dimension of the set of DAE in the indirect method. On the other hand, the solution of NLP (Nonlinear Programming) problems of great dimension or the attainment of the gradients of the objective function in the sequential method is not trivial (Feehery, 2001). The solution of OCP with inequality constraints presents an additional complexity because it demands the knowledge of the sequence and the number of constraint activations and deactivations along the trajectory. When the amount of constraints is reduced, it is usually possible to determine this sequence examining the solution of the problem without constraints. However, the presence of a large number of restrictions leads to a problem of combinatorial nature (Feehery, 2001). A particular case of great interest is the presence of a linear control variable in the Hamiltonian function. In general, no minimum optimal solution exists for such problems, unless inequality constraints in the state and/or control are specified. If the inequality constraints are linear in the control variable, it is reasonable to expect that the minimizer, if it exists, will ( ) ( ) ( ) ( ) , min , , , t f f f u t t f t o J z t t L z u t dt ? = + ? ( ) , , , 0 f z z u t = ! ( ) ( ) ( ) ( ) , , ,0o o o o z t z t u t t ? = ! , ,, , , and . (.) J (.) L ( ) . ? ? ! ( ) . f ( ) . n x ? ? ! n z z ? ! n u u ? ! always demand that the control variables are located at a point on the limits of the feasible region of control (Bryson and Ho, 1975). For this purpose, consider the following system of equations: (4) (5) with control variable given by: ( The Hamiltonian function (H) is defined as: (7) For this class of control, we have: (8) ( ) o o z t z = ( ) ( ) z F z g z u = + ! min max u u u ? ? ( ) ( ) ( ) T H F z g z u ? = + max min 0 0 0 T T T u g u g u g ? ? ? ? < ? ? = ? = ? ? > ? ? where ? is the Switching Function (Lobato, 2004). # III. # Optimization of Feed-Batch Penicillin Fermentation Process The mathematical model of the feed-batch penicillin fermentation process considered in this contribution was described and studied by San and Stephanopoulos (1989). Mathematically, this model consists of the following constraints: In this work, we formulate a robust multi-objective singular optimal control problem, based on the feedbatch penicillin fermentation process, which consists of maximizing the productivity and minimizing the operation total time, describe as: where t is the time (h), X is the biomass concentration (g/L), P is the amount of existing penicillin product (g/L), S is the substrate concentration-control variable (g/L), V is the volume of biological reactor, F is the feed rate (1666.67 L/h), µ is the growth rate and ? is the specific product formation rate. (16) ( ) 0.01 0 0 g/L dP FP X P P dt V ? = ? ? = ( ) 0 1 g/L dX FX X X dt V µ = ? = 41 g/L X ? 0.001 0.5 g/L S ? ? 0.004 1 0.0001 0.1 S S ? = + + 0.11 0.006 S S X µ = + ( ) 0 2.5E5 L dV F V dt = = ( ) ( ) max f f f P t V t t © 2020 ( )17 In order to choose a point that belongs to the Pareto Curve obtained, taking into account a multi-objective optimization strategy, the overall profit (OP) is considered. This relation is defined as (San and Stephanopoulos, 1989): (18) IV. # Multi-Objective Optimization Differential Evolution Aiming to solve the multi-objective optimization problem proposed, in this section is presented a brief review about the multi-objective optimization problem and the MODE strategy, respectively. When dealing with multi-objective optimization problems, the notion of ''optimality'' needs to be extended. The most common approach in the literature was proposed by Edge worth (1881) and later generalized by Pareto (1896). This notion is called Edge worth-Pareto optimality, or simply Pareto optimality, and refers to finding good trade-offs among all the objectives. This definition leads us to find a set of solutions that is called the Pareto optimal set, whose corresponding elements are called no dominated or no inferior. Multi-objective optimization deals with optimization problems which are formulated with some or possibly all of the objective functions in conflict with each other. Such problems can be formulated as a vector of objective functions f ( x) = [f1(x) f2(x) ? f m (x)] subject to a vector of input parameters x = [x1 x2 ? x n ], where m is the number of objectives, and n is the number of parameters. According to the criterion of Pareto, multi-objective problems have a set of trade-off solutions, where a solution may be better on objective f1 but worse on objective f2, whilst other solutions may be worse on objective f1 but better on objective f2. The literature shows a large number of multiobjective optimization techniques, although these methods have limitations when it comes to highly complex applications (Deb, 2001). Metaheuristics have established themselves as a complementary approach that can be applied even when no prior information is known about the underlying problem. The growing popularity of evolutionary algorithms in this field is mainly due to their flexibility to deal with a wide variety of multi-objective optimization problems (both numerical and combinatorial) and to their easiness of use. Also, due to their population-based nature, evolutionary algorithms can be modified such that they generate several nondominated solutions in a single run. These features have made them popular when tackling complex real world multi-objective optimization problems (Deb, 2001). In order to solve the multi-objective optimization problem, Lobato (2004) proposed the MODE strategy. This is based on the association between the Differential Evolution (DE) algorithm (Storn and Price, 1995) with two operators: ranking ordering and crowding distance. This algorithm has the following structure: an initial population of size N is generated at random. All dominated solutions are removed from the population through the operator Fast Non-Dominated Sorting. This operator calculates, for each population member, represented by x i, the number of individuals that dominate x i (generating a domination count, n i) and the set of candidates S i that are dominated by x i. Afterwards, the population is sorted into non-dominated fronts F j (sets of vectors that are non-dominated with respect to each other) as described in the following: the vectors with n i = 0 constitute the first front, F0. For every vector in the front F j (beginning with j = 0), the domination count n i of vectors of the corresponding sets S i is reduced by one. If a domination count becomes zero, the corresponding vector is put into the next nondominated front F j+1. This procedure is repeated until each vector becomes the member of a front. The remaining nondominated solutions are retained for recombination. In this step, three parents are selected at random. A child is generated from these three parents (this process continues until N children are generated). Starting from population P 1 of size 2N, neighbours are generated from each one of the individuals of the population. Those generated candidates are classified according to the dominance criterion described before and only the nondominated neighbours (P2) are put together with P 1 to form P3. The population P 3 is then classified according to the dominance criterion. If the number of individuals of the population P 3 is larger than a predefined number, the population is truncated according to the criterion defined by the Crowding Distance criterion (Deb, 2001). The crowding distance describes the density of solutions surrounding a vector. To compute the crowding distance for a set of population members the (or an arbitrary large number for practical purposes). For all other vectors, the crowding distance is calculated according to: ( V. # Effective Mean Concept Traditionally, the introduction of robustness in the multi-objective context require the consideration of new constraints and/or new objectives (relationship between the mean and the standard deviation of the vector of objective functions) and probability distribution functions for the design variables and/or objectives (Ritto et al., 2008). As an alternative to these classical formulations, Deb and Gupta (2006) extended the Effective Mean Concept (EMC), originally proposed for mono objective problems, to the multi-objective context. In this approach, no additional constraints are inserted into the original problem. Thus, the problem is rewritten as the mean of the original objectives. In this case, the robustness measure and the solution of a robust multiobjective optimization problem are defined as (Deb and Gupta, 2006): (20) (21) Where g is the inequality constraints vector and m is the number of objectives. In the present paper, the EMC is used to assess the robustness in each candidate generated by using the MODE algorithm. In this case, the original objective function vector is transformed by considering Eq. ( 21). The user needs to input the objective functions vector, the constraints vector, the design space, MODE parameters, the perturbation ? added to the vector of design variables, and the sample size N sample. # VI. # Results and Discussion In order to solve the proposed robust multiobjective singular optimal control problem, the following parameters are considered in MODE: population size (25), number of generations (200), perturbation rate (0.8), crossover rate (0.8), number of pseud-curves (10) and reduction rate (0.9). The control variable was discretized considering 5 control elements. Three cases are considered, according to the level of uncertainty: ? = 0 % (nominal solution, i.e., without uncertainty), ? = 5 % and ? = 10%. For each test case, the number of samples was equal to 50 (Nsample). Considering the parameters presented above, 25+25×200 objective function evaluations are necessary to solve the nominal case by using the MODE. In order to solve the robust cases by the MODE, 25+25×200×50 objective function evaluations are necessary. Figure 1 presents the Pareto Curve obtained by using the MODE strategy. We can observe that the increase in total operation time (tf) implies an increase in productivity. In addition, the productivity is higher for the nominal case due to higher t f values, following the robust cases. The overall profit (OP) is favored by increase of t f, as observed in Tab. 1. ( ) ( ) 1 , 1 , 1 ,max ,min 0 Where f j corresponds to the j-th objective function and m is equals to the number of objective functions. This process is executed until the total number of generations is reached. m j i j i x i j j j f f dist f f ? + ? = ? = ? ? ( ) 1 ( , ) ( ) ( ) eff y B x f x f y dy B x ? ? ? ? = ? ( ) Where x is the design variables vector, f is the objective function, f eff is the EMC applied to this function, ? is the robustness parameter, |B?|is the hyper-volume of the neighborhood in relation to the design variable x. To evaluate this integral, sample points are created randomly by using the Latin Hypercube method, in the vicinities of x. In the multi-objective context, the optimization problem is given by: As mentioned earlier, Fig. 2 presents the evaluation of OP for each individual considering nominal and robust solutions. In this case, a good diversity, in terms of individuals of the population obtained by using MODE is observed. The best individuals (see Tab. 1), in terms of the OP, are chosen to simulate the process, as observed in Figs. 3-6. In these curves, it is important to observe that, initially, the profiles are similar, due to the proximity of the feed substrate concentration of maximum value (S=0.5 g/L) to increase the cells concentration rapidly. For each value of ?, after a determine value, the feed substrate concentration reaches a value close to the minimum (S=0 g/L) to increase the product concentration rapidly. In Fig. 6 we can observe that during the first step (S?0.5 g/L), the process is not profitable due to the product concentration. # Conclusion In this paper, the MODE strategy was associated with the EMC approach to determine the feed substrate concentration in a fed-batch penicillin fermentation process. The results demonstrated that the insertion of robustness implies in the reduction of diversity of the Pareto Curve and the deterioration of the Pareto Curve in relation to nominal result. Since a systematic study introducing robustness in multi-objective optimization problems (Deb and Gupta, 2006) is not easily available, the problem studied may serve as comparison for future evaluations of other methodologies for robust multiobjective optimization. Regarding optimal robust design, the determination of robustness regions may represent a criterion for the choice of a specific point of the Pareto Curve for a possible practical implementation. However, it is important to observe that the main disadvantage of this approach is the increase of the number of objective function evaluations, which are necessary to evaluate the integral considered in the Effective Mean Concept, independently from the optimization strategy considered. Further works will be dedicated to approaches related to dynamically updating the parameters and mutation strategies of the MODE together with its parallelization to reduce the computational time. 1? (%)Operation Time Total (h)-Productivity (g/L)-overall profit ($/g) 00187.853-27613.673-1.066E+065176.652-26317.889-9.660E+0510155.605-25369.781-8.444E+05 Robust Multi-Objective Singular Optimal Control Ofpenicillin Fermentation Process * Advances in Simultaneous Strategies for Dynamic Process Optimization LTBiegler LTCervantes AMWachter Chemical Engineering Science 57 2002 * Numerical Solution of Initial Value Problems in Differential Algebraic Equations KEBrenan SLCampbell LRPetzold Classics Appl. Math. SIAM Philadelphia 1996 * Applied Optimal Control AEBryson YCHo 1975 Hemisphere Publishing Washington * Singular Control in Bioreactors ACCosta 1996 Brazil PEQ/ COPPE/ UFRJ, Rio de Janeiro M.Sc. Thesis in Portuguese * Multi-Objective Optimization using Evolutionary Algorithms KDeb 2001 John Wiley & Sons Chichester, UK * Introducing Robustness in Multiobjective Optimization KDeb HGupta Evolutionary Computation 14 4 2006 * FYEdgeworth 1881. Mathematical Physics, P. Keagan * Dynamic Optimization with Path Constraints WFFeehery 2001 MIT PhD thesis * Non-singular Optimal Control for Fed-Batch Fermentation Processes with a Differential-Algebraic System Model PCFu JPBarford Journal Process Control 3 1993 * Optimal Substrate Feeding Policy for a Fed Batch Fermentation with Substrate and Product Inhibition Kinetics JHong Biotechnology and Bioengineering 28 1986 * Multi-objective Optimization to Engineering System Design FSLobato 2008 Uberlândia-MG, Brasil Tese de Doutorado, Universidade Federal de Uberlândia * General Characteristics of Optimal Feed Rate Profiles for Varies Fed-Batch Fermentation Processes JMModak HCLima YJTayeb Biotechnology and Bioengineering 28 1986 * Simple Non-singular Control Approach to Fed-Batch Fermentation Optimisation JMModak HLim Biotechnology and Bioengineering 33 1989 * V1896Pareto ;Cours D'economie Politique II FRouge Lausanne * Timoshenko Beam with Uncertainty on the Boundary Conditions TGRitto RSampaio ECataldo Journal of the Brazilian Society of Mechanical Science and Engineering XXX 4 2008 * Optimization of Fed-Batch Penicillin Fermentation: A Case of Singular Optimal Control with State Constraints K.-YSan GStephanopoulos Biotechnology and Bioengineering 34 1989 * Differential Evolution -A Simple Evolution Strategy for Fast Optimization RStorn KPrice Dr. Dobb's Journal 22 4 1995 * Modelling and Optimal Control of Fed-Batch Processes using a Novel Control Affine Feed Forward Neural Network ZXiong JZhang Proceedings of the 2002 American Control Conference the 2002 American Control ConferenceAnchorage, AK, USA 2003