Methods and Tools for Robust Optimal Control of Batch Chromatographic Separation Processes

This contribution concerns the development of generic methods and tools for robust optimal control of high-pressure liquid chromatographic separation processes. The proposed methodology exploits a deterministic robust formulation, that employs a linearization of the uncertainty set, based on Lyapunov differential equations to generate optimal elution trajectories in the presence of uncertainty. Computational tractability is obtained by casting the robust counterpart problem in the framework of bilevel optimal control where the upper level concerns forward simulation of the Lyapunov differential equation, and the nominal open-loop optimal control problem augmented with the robustified target component purity inequality constraint margin is considered in the lower level. The lower-level open-loop optimal control problem, constrained by spatially discretized partial differential equations, is transcribed into a finite dimensional nonlinear program using direct collocation, which is then solved by a primal-dual interior point method. The advantages of the robustification strategy are highlighted through the solution of a challenging ternary complex mixture separation problem for a hydrophobic interaction chromatography system. The study shows that penalizing the changes in the zero-order hold control gives optimal solutions with low sensitivity to uncertainty. A key result is that the robustified general elution trajectories outperformed the conventional linear trajectories both in terms of recovery yield and robustness. Processes 2015, 3 569


Introduction
Isolation of a high-purity target component from a multicomponent mixture is of significant importance in the pharmaceutical and biochemical industries [1].In the clinical or commercial-scale manufacturing of human therapeutic proteins, high-pressure liquid chromatography (HPLC) is an essential process operation to achieve the high purity requirements for biopharmaceutical drugs [2,3].Due to inflexibility in the conventional regulatory paradigm, applying process operational changes to an approved process design in the biopharmaceutical community has been historically prohibited without filing for a new approval [1].The United States Food and Drug Administration (US FDA) has published a series of new guidelines to promote flexibility for innovation while appropriately managing variability in the process performance attributes for continual improvement via process analytical technology (PAT) and quality-by-design (QbD) principles [4,5].QbD promotes improved process understanding during process and product development and building quality in the design instead of testing for quality.This can be achieved via correlative, causal, or mechanistic knowledge and at the highest level via first principle models.Hence, PAT has been gaining a lot of momentum in the biopharmaceutical community due to the potential for continuous real time quality assurance resulting in improved operational control and compliance [6].
In the QbD approach, critical quality attributes (CQA) of pharmaceutical products are defined that assure desired clinical performance, and then a manufacturing process is designed to consistently meet these product attributes, thus assuring product quality [7].Process characterization is conducted to identify the impact of process parameters on the products CQA's [8], which is then used to define a process design space.A design space is defined as "the multidimensional combination and interaction of input variables and process parameters that have been demonstrated to provide an assurance of product quality" [9].Acceptable ranges for process parameters and input variables are documented in the regulatory filing, and working within these ranges is not deemed to be a change from normal operating conditions.The expected benefit of the QbD approach is an increase in the assurance of product quality, and in turn, the FDA will allow manufacturers greater flexibility to operate with lower regulatory burden, enabling continuous process improvement, as well as greater robustness [10].
The current approach to process characterization within the QbD framework is concerned with validating that the defined design space can cope with the process variability experienced during normal operation [10].In this context, the FDA guidance encourages the application of mechanistic models to improve process understanding, based on fundamental knowledge of the underlying physico-chemical phenomena correlating the nonlinear interdependence of process parameters and inputs to the resulting value of the CQA's.Model-based approaches to sensitivity and robustness analysis have been successfully applied to a diversity of chromatographic separations, see e.g., Mota et al. [11], Degerman et al. [12], Westerberg et al. [13], Borg et al. [14] and the references cited therein.The aforementioned robustness assessment investigations have exclusively been founded on stochastic methodologies, in which inherent process parameters and inputs have been randomly generated from a uniform or a truncated normal distribution, to determine the magnitude of variability of the CQA's.The result is a probability density distribution of the output, which can be used to determine the risk of batch failure and the process capability as well as to identify inputs that can be used to design a control strategy.In this context, the main objective of this study is to exploit recent advances in deterministic robustness strategies in conjunction with open-loop optimal control.

Solvent Composition Trajectory Elution Strategies
Solvent composition trajectory elution is widely applied in both analytical and preparative chromatography, and refers to a continuous change in the mobile phase during separation [15].By this means, improved resolution of highly complex sample mixtures can be obtained in a much shorter time than could be expected in isocratic elution conditions.Reproducible selectivity in chromatographic separation of closely related product impurities requires optimized elution mode and fractionation interval endpoints [16].The current state-of-the-art methodology for optimization of HPLC separations has been limited to linear [17,18], concave/convex [19] and step elution trajectories [20].A novel open-loop optimal control strategy for simultaneous optimization of general elution trajectories and target component fractionating decisions was recently developed by Holmqvist and Magnusson [21], Holmqvist et al. [22].The present investigation is a direct continuation of these studies and the overall objective is to extend the open-loop optimal control strategy with a robust worst-case formulation and benchmark the robustness of general elution trajectories, with respect to target component purity, with that of the conventional linear trajectories.
In the context of the QbD framework, this implies relating the impact of a critical process parameter (CPP) to the difference between the quality requirement and the CQA at the nominal operating point, i.e., assess the level of reduction in the robustness safety margin.A CPP is formally defined as "a process parameter whose variability has an impact on a critical quality attribute and, therefore, should be monitored or controlled to ensure the process produces the desired quality" [23].A quantitative value of the sensitivity of the CQA with respect to CPPs can be obtained by scaling the discrepancy between the CQA value at the nominal operating point and the requirement by perturbation of the CPP values of 1.0 standard deviation in the direction that gives the worst-case value of the QCA [13]: The higher the value, the more likely the parameter is to cause batch failure.Thus, it is noteworthy that the open-loop controlled system targeting the captured amount of the target component with respect to the fractionating interval endpoints yields that the inequality constraint imposed on target component purity is active, i.e., the resulting purity equals its lower requirement [16].Consequently, there is no inherent safety margin incorporated in the open-loop optimal control problem and to surmount this shortcoming calls for the formulation of a robust counterpart problem [24].

PDE-Constrained Dynamic Optimization
In this study, the realistic multi-component system dynamics required for analysis were generated by numerical solution of the reaction-dispersive model [25].This model is governed by a set of mass-balance partial differential equations (PDEs), with a modified Langmuir isotherm and experimentally validated kinetics.In this regard, the advantages of this robust optimal control framework are highlighted through the solution of a challenging ternary mixture separation problem, with human insulin (insulin aspart, desB30 insulin and insulin methyl ester) and the intermediately eluting component as the target, for a hydrophobic interaction chromatography (HIC) [26,27] system.
The proposed model-based methodology for robust optimal control implies formulating and solving a large-scale dynamic optimization problem (DOP) constrained by PDEs [28].Optimal design, optimal control, and parameter estimation of systems governed by PDE give rise to a class of problems known as PDE-constrained optimization [29,30].The size and complexity of the discretized PDEs often pose significant challenges for contemporary optimization methods.There exists two generic direct methods to transform the infinite DOP into a finite dimensional nonlinear program (NLP); sequential and simultaneous [31,32].Optimization studies of HPLC separations carried out in the literature have exclusively been using sequential methods, where the system dynamics constraint is handled by embedded numerical integrators, and where the PDEs are approximated using the method-of-lines [33,34] and Galerkin finite element or finite volume methods.Both gradient-based NLP solvers (e.g., sequential quadratic programming and interior-point method) [35][36][37] and those based on trajectory free heuristic approaches (e.g., generic algorithms and simulated annealing) [18,38,39] have been used successfully.The proposed methodology in this paper is based on a developed simultaneous method where both the control and state variables are fully discretized in the temporal domain using orthogonal collocations on finite elements [31].In order to reduce the size of the resulting NLP, the PDE system was approximated using an adaptive, high order finite volume weighted essentially non-oscillatory (WENO) scheme [40][41][42].The NLP is subsequently solved using the primal-dual interior point method IPOPT [43] and algorithmic differentiation (AD) techniques [44].
Jacobians are required by gradient-based NLP algorithms and simulation.Access to accurate Jacobians often improves the performance and robustness of algorithms, and in addition, efficient implementation of Jacobian computations can reduce the over-all execution time.In the context of HPLC system modelling, Püttmann et al. [45] compared the accuracy of Jacobians, with respect to an intrinsic model parameter, computed using first-order finite differences (FD) and AD.It was shown in that study the AD approach outperformed the FD approach in terms of both accuracy and over-all execution time for realistically sized chromatography models when computing forward sensitivities of the DAE system, i.e., the linearization of the original system with respect to an intrinsic model parameter.The AD approach in [45] is extended in this paper to include computations of Hessians required by the NLP solver but also addresses the construction of metrics required in the robust counterpart formulation.

Nonlinear Robust Optimization of Uncertain Dynamic Systems
As outlined above, uncertainty propagation of CPPs through nonlinear model equations is reflected in terms of QCAs (e.g., constraints and objectives).Methodologies developed to cope with the complexity of optimization problems under uncertainty can be categorized into stochastic programming, probabilistic (chance-constraint) programming and fuzzy programming [46].
The stochastic programming formulations assume that the probability distributions governing the uncertain parameters are available or can be estimated from the existing data.Here, reliability requires feasibility for all the outcomes of uncertain parameters, whereas chance-constrained programming requires feasibility of solutions with at least some probability specified on constraints having uncertain parameters.The main advantage of the chance-constrained programming technique is the emergence of a deterministic equivalent problem [47].In the deterministic formulations, the nature of the uncertainty is assumed to be bounded to a prespecified set, and reliability requires that all constraints are satisfied in all possible worst case situations.In any case, including more robustness gives rise to a lower performance, and hence, robustness can be interpreted as an additional and conflicting objective [48].
In this study we concentrate on the case that the dynamic system governed by PDEs is affected by a time variant input with uncertainties.To efficiently and accurately generate optimal elution trajectories in the presence of uncertainty we consider a deterministic robust formulation, that employs a linearization of the uncertainty set, proposed in the article series by Logist et al. [49], Houska and Diehl [50], Houska et al. [51].The main objective is to robustly regard inequality state constraints.For this aim, the proposed methodology utilizes Lyapunov differential equations to compute variance-covariance matrix functions and facilitates the tractable robust counterpart formulation for optimal control problems.

Aim and Scope
External process disturbances are inevitably present and the robustness of the open-loop controlled system with respect to CPP uncertainties, e.g., guaranteeing that QCA requirements are not violated, is of the highest importance in biopharmaceutical process industry.It is, however, noteworthy that there is a drastic increase in computational complexity while moving from the nominal problem to its robust counterpart.The main contribution of this paper is therefore to establish a methodology and a computationally efficient framework for robust optimal control yielding robust general elution trajectories, with respect to target component purity, for batch chromatographic separation processes.The study presented here had three main objectives: i) To formulate a nominal open-loop optimal control problem and its robust counterpart, including penalty on the difference on the zero-order hold control discretization, to find optimal control trajectories that remain feasible for all perturbations from the uncertainty set.ii) To develop a numerically efficient framework for simulation of a Lyapunov differential equation, enabling the computation of the robustified inequality constraint margin in the counterpart problem.iii) To assess the robustness of general elution trajectories, for both batch elution and cyclic-steady-state operation, and to benchmark with that of the conventional liner trajectories.

Outline of the Paper
The remainder of this paper is structured as follows: Section 2 presents the HPLC separation system.Section 3 outlines the process model and the spatial discretization scheme.Section 4 describes the open-loop and robust counterpart problem formulation.Section 5 outlines the methods and tools for dynamic simulation and optimization.Section 6 presents the primary results, and Section 7 contains concluding remarks and perspectives for future research.

Process Description
The HPLC separation system operated in batch mode is characterized by subsequent pulse injections of the sample mixture and relies on pumps to pass a pressurized liquid solvent containing the sample mixture through a column packed with a solid adsorbent material [52].A simplified P&ID of a standard HPLC system under consideration is illustrated in Figure 1.The system facilitates implementation of zero-order hold elution trajectories, u, through control of the volume fraction of buffers with different solvent strengths [53].The resulting time-variant mobile phase modifier concentration, c mix,S , is obtained in the mixing unit and enters the separation column with the system volumetric flow rate Q.Moreover, the feed, with concentration c load,α and α ∈ {A, B,C}, is injected as a rectangular pulse with duration ∆t load via a multiport valve.The eluting concentration trajectory is analyzed with a UV detector at the column outlet [54], and a fractionating valve separates the mixture into its components [36].The fractionating interval manipulating variables [τ 0 , ∆τ] govern the maximum recovery yield while still fulfilling the constraint imposed on purity of the target component fraction [16].The time-invariant operating parameters [ Q, c load,α , ∆t load ] are summarized in Table 1 and the optimal control problem optimization variables include the piecewise constant controls of the parameterized elution trajectory, u i and ∀i ∈ [1..N u ], and [τ 0 , ∆τ]. (4) Simplified P&ID of the HPLC system.Elementary system component description: (1) buffer mixing unit, (2) high-pressure switching valve unit, (3) HPLC separation column, (4) UV detector, and (5) fractionating valve.
a Pe is assumed to be constant and equal to 0.50, i.e., molecular diffusion is neglected [55].

Control Signal Uncertainties
The control signal uncertainty is assumed to arise from off-set volumetric fractions in the buffer mixing and initial off-set in the individual buffer concentrations.Robustness is often associated with the ability for a feedback controller to compensate for disturbances and uncertainties.However, chromatographic separations cannot easily be controlled by conventional control strategies due to their complex dynamics with extremely long time delays, spatially distributed properties, and switchings [36].For these reasons, the open-loop optimal control problem, which yields that the inequality constraint imposed on target component purity is active, is extended with its robust counterpart to cope with time-variant variability in u.It is noteworthy that u has a corresponding actuator system, i.e., the control signal is in fact a set-point to the actuator system which includes a low-level controller.The limits for the low-level controller is approximated by absolute and rate limits on the change in the zero-order hold control, ∆u.It is expected that a bang-bang solution would not be robust in the presence of time-variant control signal uncertainty, since even small input perturbations would highly influence the eluting concentration trajectories, with degraded purity of the intermediate eluting target component as a result.In order to cope with bang-bang elution trajectories, a quadratic cost on the difference of the piecewise constant control flows, ∆u, is incorporated in the optimal control problem which will influence the smoothness of u and consequently enhance the robustness, however, to the expense of the recovery yield.

Mathematical Modeling
The governing one-dimensional equations in the reaction-dispersive model [25] of the mobile and stationary phase-Derived under the assumptions of infinitely fast diffusion into the particles and rate-limiting adsorption kinetics-Defined in the spatial, z ∈ [z 0 , z f ], and temporal, t ∈ [t 0 ,t f ], domains are: where c α and q α are the mobile and stationary phase concentration of component α ∈ {A, B,C, S}, v int denotes the interstitial velocity of the fluid, D app the apparent dispersion coefficient, and ε c and ε p the column and particle void fractions.In this study, the self-association isotherm [56] was used to describe the hydrophobic interaction in Equation ( 2).Here, c S denotes the concentration of the non-absorbing modifier (i.e., ∂ q S /∂t : = 0), k kin the kinetic rate constant, H 0,α the Henry's constant, γ α the solvophobicity parameter, ν the binding charge ratio, q max the maximum concentration of adsorbed components, and K eq the equilibrium constant.Equation ( 1) is complemented with Danckwerts boundary conditions: where c load,α is the injected load concentration and Π(t,t 0 , ∆t load ) ∈ {0, 1} a rectangular function which is 1 in the temporal horizon [t 0 ,t 0 + ∆t load ].The dynamics of the modifier concentration in the upstream mixing tank, c mix,S , is governed by: dc mix,S dt = 1 where u is the set-point elution trajectory and τ mix the residence time.It is noteworthy that the control u is the absolute modifier concentration and not the relative volumetric fraction of buffers with different solvent strengths as previously outlined in Section 2.
The set of parameters associated with the adsorption's dependency on mobile phase concentrations and the properties of the solvent and the stationary phase was calibrated to lab-scale experimental data, to reproduce the behavior in the studied system by means of the inverse method [57][58][59].A detailed description of the experimental design and the materials used is outlined in [60], and the least-square estimates of the adsorption isotherm kinetics are listed in Table 2.

Table 2. Kinetic parameters of the Langmuir self-association adsorption model Equation (2).
The feed composition, c load,α , is given in fractional weights for each insulin component and the total feed concentration is 2.0 × 10 1 (moL

Spatial Discretization
The equations of the spatially distributed HPLC column model that describe the mobile and stationary phase state dynamics (Equations ( 1) and ( 2)), constitute a system of non-linear PDEs.In this study, the PDE system was approximated using the method-of-lines [33,34] and high-order finite volume WENO scheme [40][41][42] on a uniform mesh where z j = j∆z is the discretized spatial coordinate and j ∈ The main advantage of this scheme is its capability to achieve high-order formal accuracy in smooth regions while maintaining stable, nonoscillatory, and sharp discontinuity transitions.In contrast to other high-order spatial discretization schemes that make use of flux limiters [61], such as the MUSCL scheme [62] and the TVD scheme [63], the WENO scheme does not require tuning of parameters and more importantly the discretization scheme is twice continuously differentiable; a necessity when employing Newton-based methods for solving NLPs.Although a high-order WENO scheme may use several times more CPU time than the aforementioned schemes [64], which are usually second-order accurate in the smooth part of the solution, it is still computationally advantageous for the purpose of this study where the DOP is transcribed into an NLP using direct collocation in time.Hence, less grid points are required to spatially resolve complicated smooth structures with small numerical dissipation, and thereby reducing the size of the resulting NLP.

Problem Formulation
The cyclic operation of HPLC separation processes in batch elution mode imposes adequate retention time repeatability by conditioning the column to the initial modifier concentration prior to the subsequent injection [65,66].In the scope of chromatographic separations, cyclic-steady-state (CSS) prediction has previously only been applied to continuous simulated moving bed (SMB) [11,18,67] and multicolumn counter-current solvent gradient purification (MCSGP) [68] processes.Hence, the CSS operation implies computation of a limit cycle dynamic solution, subject to initial conditions for the state variables, that conform to the periodic boundary condition over the processing cycle.It is however noteworthy that in order for the comprehensive CSS description to be experimentally feasible, a model representation of the most retained component in the loaded sample mixture is required.This is vital for prohibiting potential stationary phase accumulation of otherwise irreversibly adsorbing components, resulting in a gradual decrease in column adsorption capacity.For these reasons, the robustification methodology considered in this study was applied to two inherently different optimal control strategies: i) Optimal control strategy I is exclusively concerned with the batch elution operation mode, where the limit cycle criteria are relaxed and the upper fractionating interval endpoint bounds the temporal horizon, i.e., t f = τ 0 + ∆τ.ii) Optimal control strategy II is concerned with the comprehensive CSS operation mode, where the temporal horizon includes the column regeneration and re-equilibration modes.
Hence, optimal control strategy I that considerably reduces the DOP complexity, is preferable when the complete sample mixture composition is unknown and is extensively applied in open-loop optimal control studies on non-isocratic batch elution chromatography, see, for example, Damtew et al. [19], Osberghaus et al. [69] and the references cited therein.

Cyclic-Steady-State Criteria Formulation
The computation of CSS solutions over the temporal horizon t ∈ [t 0 ,t f ] requires that Equations ( 1) and ( 2) are augmented with additional criteria governing that the state at the initial time is retained at the end of the cycle.Accordingly, the limit cycle criteria can be classified into terminal equality constraints: and initial conditions satisfying: where δ load,α = c load,α v int A c ∆t load is the total injected sample amount.Specifically, Equation (6a) states the flux of the most retained component α = C at the column outlet, z = z f , equals the total loaded sample amount, δ load,α , at t = t f .This criterion, thereby, ensures that all components are completely eluted at the end of time horizon.The terminal equality constraints Equation (6b)-(6d) govern that the modifier concentration, c S , at every column position z ∈ [z 0 , z f ] as well as the concentration in the mixing unit, c mix,S , are consistent at the initial and terminal times.Moreover, Equation (6e) and (6f) are supplementary initial conditions, enforcing the initial modifier concentration at t 0 .Equation (6a)-(6d) are introduced in the DOP whereas Equation (6e) and (6f) are introduced as initial values in the governing mobile phase transport equations.

Open-loop Optimal Control Problem Formulation
In chromatographic separations, there are several incommensurable objectives which require a trade-off to ensure satisfactory design, see e.g., Guiochon et al. [70] for a thorough review.
Typical objective functions are the recovery yield, Y α , the production rate, P α , and the recovery concentration, c α , with respect to the target component collected in the fractionating interval.In any case, the elution trajectories at the column outlet, c α (t, z f ) for α ∈ {A, B,C}, form the basis for evaluating the incommensurable objective functions.These are defined here as: and subject to initial values Y α (t 0 ) = 0, P α (t 0 ) = 0, and c α (t 0 ) = 0.Moreover, Π(t, τ 0 , ∆τ) is a smooth rectangular function in the fractionization horizon [τ 0 , τ 0 + ∆τ].It is noteworthy that both P α , defined as the amount of target component collected per cycle time scaled to the size of the column, and c α , defined as the integral mean value of the target component concentration over the fractionating horizon, are intrinsically governed by Y α .A weighted sum scalarization method is frequently used to combine the objectives in Equations ( 7)-( 9) into a single performance index.The procedures' deterministic nature allows the use of direct NLP routines and ensures an accurate and efficient Pareto set generation [22].
For the purpose of this study, we solely consider the yield of the intermediately eluting component α = B as target in the optimal control problem.Hence, the aim is to derive an optimal elution trajectory, u, and fractionating interval endpoints, [τ 0 , τ 0 + ∆τ], that maximize Y α (t f ) while fulfilling the constraint imposed on purity of the target component fractionization: where the numerator of the right hand side represents the captured amount of the target component in [τ 0 , τ 0 + ∆τ] and the denominator represents the total amount captured.Hence, X α (t f ) is incorporated in the DOP as a terminal inequality constraint, X α,L −X α (t f ) ≤ 0, with an assigned lower purity requirement X α,L .Given the optimization specifications in Equations ( 7) and ( 10) the open-loop optimal control problem, with differential algebraic constraints [31,71], over the temporal domain [t 0 ,t f ] may now be formulated as: In addition to the free operating parameters, i.e., the discrete control signal u and the free time-invariant parameters (τ 0 , ∆τ), the optimization variables include the state variables x(t) = (c α (t, z j ), c S (t, z j ), c mix,S (t), q α (t, z j ),Y α (t)) for α ∈ {A, B,C} and the algebraic variable X B .F denotes the implicit DAE system resulting from the spatially discretized PDAE system.Moreover, g e denotes the terminal equality constraints and assembles the limit cycle criteria defined in Equation ( 6).Time-invariant bounds on the fractionation variables are introduced in Equation (11f) and absolute and rate limitations of the control are enforced in Equation (11e).As outlined in Section 2.1, Equation (11a) is augmented with a quadratic penalty on the differences of the piecewise constant controls, ∆u, discretized with N u segments in the temporal domain [t 0 ,t f ].This will influence the smoothness of u, and the weight R is used as an optimization variable in the robust counterpart problem in order to find optimal robustified elution trajectories.

Robust Counterpart Problem Formulation
The solution of Equation (11) will yield optimal open-loop trajectories that lie on the boundary of the feasible region.This means that even the smallest of unmodeled disturbances can cause infeasbility, leading to failed batches.To account for this, we discuss how a robust counterpart of Equation ( 11) can be formulated and then approximated based on linearization to reduce computational complexity.The resulting formulation is still not computationally tractable, but we present how it nevertheless can be used to find robust solutions using bilevel optimization.
In parts of this section we will assume that the system dynamics in Equation (11b) are described by an explicit ordinary differential equation (ODE) system rather than a DAE.Since the considered DAE is of index one, the transformation to an ODE is straightforward and is handled by the framework discussed in Section 5.1.2.The uncertainties discussed in Section 2.1 are modeled as a disturbance w that is added to the control variable u, whose L 2 norm is bounded by Γ.The constraint that we want to make robust with respect to these uncertainties is the purity constraint X B (t f ) ≥ X B,L .The considered robust counterpart of Equation ( 11) is thus essentially to modify X B (t f ) in the purity constraint to instead be the minimum of X B (t f ) over all disturbances w with L 2 norm less than Γ (A rigorous formulation of the considered robust counterpart involves transforming the problem to Mayer form and eliminating the state variables by utilizing the existence of a unique solution of the ODE.The result is verbose and neither enlightening nor useful for computations, and is thus omitted).The result is a problem that is not computationally tractable due to the nonlinear dynamics, even if the state dimension would have been small.To reduce the computational complexity, we utilize the linearization approach presented in [49].We first solve the original problem Equation (11) to obtain nominal trajectories x * and u * .We then introduce the deviation δx caused by the disturbance w, leading to the perturbed trajectory x = x * + δx.The deviation δx can be approximately determined by linearizing the dynamics around the nominal trajectories x * and u * .This leads to the uncertain linear time-varying system for determining the deviation corresponding to w, where The utilization of the above linearization, which is valid for small disturbances, and [Theorem 3.1] [72] yields that the robust counterpart can be equivalently formulated by the introduction of the back-off term Γ C(t)P(t)C T (t), where and P : [t 0 ,t f ] → S n x + satisfies the Lyapunov differential equation The back-off term is used to augment the purity constraint Equation (11c) to get This back-off term can be regarded as an additional safety margin required for guaranteeing the meeting of the constraint.However, to compute it, we need to introduce P, givings us n x (n x + 1)/2 extra states.The result is the following approximate robust counterpart of Equation (11).
The discussed robustification also has a stochastic interpretation [72].Rather than considering disturbances with bounded norm Γ, we can instead consider stochastic disturbances.If w is a Gaussian white noise process with an identity covariance matrix, then P is a linear approximation of the covariance matrix of the state x.Consequently, the covariance matrix associated with the output function h is approximated by C(t)P(t)C T (t).While the discussed approximate robust approach is computationally tractable for medium-scale ODE-constrained DOPs, it is not for our original PDE-constrained DOP.For these reasons, the robust counterpart problem is cast in the frame of bilevel optimal control where the upper level concerns forward simulation of the Lyapunov differential equation, and the nominal NLP augmented with the robustified inequality constraint margin is considered in the lower level: w.r.t.
x : It is noteworthy that while the nominal optimal control problem Equation ( 11) is feasible, the robust counterpart formulation Equation (19) might not be feasible for large disturbances as the constraints are required to be satisfied.The upper-level static optimization problem (Equations (19a)-(19d)) with respect to R is solved using a Nelder-Mead Simplex algorithm, whereas the lower-level DOP (Equations (19e)-(19j) constrained by the DAE system dynamics is transcribed into an NLP using direct collocation, as described in Section 5.2.This strategy avoids the full simultaneous discretization of the Lyapunov states, which are too numerous to be treated with direct collocation.Finally, it is important to clarify that the Lyapunov state simulation and the subsequent NLP execution is repeated in the lower level until an optimal solution is obtained for a fixed R, i.e., Equation (19g) is satisfied to a small specified tolerance.

Methods and Tools
This section outlines the languages and tools used to generate the simulation and optimization results.We first present the general software framework and then discuss the methods in the framework that allow for the efficient solution of Equation (19).

The JModelica.org Toolchain
JModelica.org[73] is an open-source platform for simulation and optimization of models described in either the Modelica [74] or Optimica language [75] or as Functional Mock-up Units (FMUs) following the Functional Mock-up Interface (FMI) [76].The platform consists of a Modelica and Optimica compiler meaning that the input is a model described in either of the two languages and the output of the compiler is either an FMU or a symbolic representation of an optimization problem.Modelica is an equation based, high-level language for describing complex physical systems.Optimica is an extension which additionally allows for specifying dynamic optimization problem with a cost function and constraints.

Optimization
The dynamic optimization framework in JModelica.org is centered around CasADi [77] and an overview is shown in Figure 2. CasADi (Computer algebra system with Automatic Differentaion) is an open-source, low-level tool for efficiently computing derivatives using AD and is tailored for dynamic optimization.The user encodes the dynamic optimization problem using Modelica and Optimica.The JModelica.orgcompiler then transfers a symbolic representation of the problem to CasADi Interface [78], after performing symbolic transformations such as index reduction.CasADi Interface serves as a three-way, symbolic Python interface between the optimization problem, user, and numerical algorithm.The main numerical algorithm for dynamic optimization in JModelica.org is based on direct, local collocation and is described in Section 5.2.The algorithm is implemented using CasADi to transcribe the DOP into an NLP.First-and second-order derivatives of the NLP functions can then be conveniently and efficiently computed using CasADi's algorithmic differentiation while preserving sparsity.The NLP is then finally solved by the open-source, interior-point method IPOPT [43].

Simulation
The simulation capabilities in the JModelica.orgplatform is based on FMUs and collected in the standalone Python package PyFMI.It is a package for interacting with models following the FMI and designed to provide a high-level, easy to use environment for working with FMUs.A key feature is the connection to Assimulo [79] which provides capabilities for performing simulations using state-of-the art solvers such as CVode from the Sundials suite [80].PyFMI together with Assimulo provides an efficient and flexible environment for simulation of FMUs.An overview of the simulation capabilities is shown in Figure 3.The FMUs are typically created from a Modelica model using the JModelica.orgcompiler.However, as FMI is an open standard, any tool that support generation of FMUs can be used.

Transcription of the Dynamic Optimization Problem
The optimization problem Equation ( 11) falls within the class of dynamic optimization problems of the general form: min.
where y is the algebraic variable and p is the free time-invariant parameters.The objective Equation (20a) is a typical optimal control Lagrange term together with a quadratic penalty on ∆u i , where R ∈ R n u ×n u acts as a weight and is typically diagonal.The optimization variables are the free operating parameters-The discrete control signal u i and free time-invariant parameters p-And the trajectories x and y.The trajectories are determined by the free operating parameters via the implicit DAE system in Equation (20b).
In this section we present an overview of the method employed by JModelica.orgto solve Equation (20) and in particular the low-level problem in Equation (19).A more detailed description is available in [81].The method is based on direct local collocation [31,32].The fundamental idea is to discretize the differential equations, thus transforming, or transcribing, the dynamic optimization problem into a finite-dimensional NLP, which then is solved by IPOPT.The discretization scheme is based on collocation methods, which correspond to special cases of implicit Runge-Kutta methods and are also commonly used for numerical solution of DAE and stiff ODE systems [82].
We start by dividing the optimization time horizon into n e elements.Within each element i ∈ [1.
.n e ] the trajectories x and y are approximated using low-order polynomials, which are called the collocation polynomials for that element.The collocation polynomials are formed by first choosing n c collocation points, which are chosen to be the same for all elements.Let t i,k be collocation point k ∈ [1..n c ] in element i and (x i,k , y i,k , u i,k ) denote the value of (x(t i,k ), y(t i,k ), u(t i,k )).The collocation polynomials in element i are then constructed by interpolating the values (x i,k , y i,k ).The polynomials for ẋ are obtained by differentiating the polynomials for x.The control variable values u i,k are obtained by constructing a map σ : [1..n e ] × [1..n c ] → [1..N u ] such that u σ(i,k) = u(t i,k ).There are different schemes for choosing the collocation points τ k , with different numerical properties, in particular regarding stability and order of convergence.We use Radau collocation, which always places a collocation point at the end of each element, and the rest are chosen in a manner that maximizes convergence order.
As decision variables in the transcribed NLP we choose the system variable values in all the collocation points, ẋi,k , x i,k , y i,k , the discrete control signal, u i , the state at the start of each element, x i,0 , and the free parameters p.The transcription of Equation ( 20) results in the NLP min.
The transcription of the Lagrange term in Equation (20a) utilizes Gauss-Radau quadrature within each element to approximate the integral by a sum using the quadrature weights ω k and element lengths h i .Equation (20b) and (20c) are transcribed into Equation (21b) and (21c), respectively, by enforcing them only in each of the collocation points, rather than during the entire time horizon.The constraints Equation (20d) and (20e) are finite in number and are thus straightforward to transcribe into Equation (21d) and (21e), respectively, using the collocation point values.
Since the states need to be continuous (but not differentiable) with respect to time, the new continuity constraint Equation (21f) needs to be introduced.Because we use Radau collocation, where no collocation point exists at the start of each element, this also requires the introduction of the new variables x i,0 , which represent the value of the state at the start of element i.Finally, we introduce Equation (21g) to capture the dependency between x and ẋ, which is implicit in Equation (20).The weights α l,k are obtained from differentiating the collocation polyonimals for the state and related to the butcher tableau of the Runge-Kutta method that corresponds to the collocation method.

Simulation Setup and Coupling of the Lyapunov Equation
In this study, the process model is implemented in Modelica and compiled into an FMU, for simulation, using the tool Dymola [83].Any tool that is able to transform the Modelica model into an FMU can be used, provided that support for directional derivatives is included.Hence, the process model is transformed and exposed as a system of ordinary differential equations, ẋ(t) = f (t, x(t), u(t), τ 0 , ∆τ), This transformation step is performed by the exporting tool.As the process model is stiff an implicit solver is needed to simulate the problem which in turn requires the Jacobian matrix of Equation (22).The Jacobian can either be approximated using finite differences, symbolic differentiation or using algorithmic differentiation.However, in order to simulate this model efficiently and accurately, access to the Jacobian, either symbolically or via algorithmic differentiation is crucial.In the FMU, directional derivatives are available and using these, the partial derivatives can be constructed (i.e., the Jacobian).The directional derivatives are defined as, where g(z) is a set of equations, in our case, either Equations ( 22) or (23), z is either the states x or the inputs u and v is the seed vector.From the directional derivatives, the partial derivatives can directly be computed by a sequence of calls with v replaced by unit vectors.However, as discussed below, there are more efficient approaches.For efficiency, FMI provides structural information in the sense that the dependency information between the functions f i and x j , u k are available.Utilizing this information, multiple columns can be constructed for each call to the directional derivative method due to that all equations are not dependent on all of the states.Finding the minimal set of calls to the directional derivative method for a general problem is difficult.However, using the algorithm by Curtis-Powell-Reid [84] one can find a set that substantially reduces the amount of calls in a straight-forward manner.In the considered case, the number of calls was reduced by 90% by applying the algorithm.Additionally, the structural information directly corresponds to the sparsity pattern which is shown in Figure 4.The pattern indicates that the model is sparse and using a sparse solver for the problem is thus beneficial.Using both the structural information and the directional derivatives, the process model, as an FMU, can efficiently and accurately be simulated.However, our interest is not solely in the direct process behavior, it is also in the matrix valued Lyapunov states, P(t), together with A(t) and B(t) defined in Equation ( 14).These equations significantly increase the problem dimension making it crucial that both the sparsity is used during the simulation and that the structure is used during the assembly of the partial derivatives.
As the solver does not support matrix differential equations, Equation ( 16) needs to be re-written into vector form.Using the notation vec(A(t)) for vectorization of A(t) by stacking the columns into a single column vector.Then, the Lyapunov equation can be written as, where ⊗ is the Kronecker product.Both the Lyapunov equations and the process equations needs to be solved together.This has been realized by extending PyFMI to support adding additional equations to the original problem of simulating an FMU.Additionally, Assimulo has been extended to be able to handle sparse representation of the Jacobian together with a connection to the sparse solver in CVode.
In Figure 3, the connection between the packages, the solver and the Lyapunov equations are shown.Putting it all together, the original equations appended with the Lyapunov equations, can be written as, together with the Jacobian, In Figure 5 the sparsity pattern of the coupled system is shown.Due to the substantially increased problem size it is even more important that A(t) can be computed efficiently and the full Jacobian represented as a sparse matrix.Finally, the quantity C(t)P(t)C T (t) governing the back-off term in Equation ( 17) is also computed.

Results and Discussion
In this section we discuss the solution of the nominal open-loop optimal control problem Equation ( 11) and its robust counterpart problem Equation (19).The uncertainty aspect is introduced via the zero-order hold elution trajectory which has a direct effect on the target component purity inequality constraint Equation (10), and hence, a robustification of this constraint is focused on.It is noteworthy that variability in CQAs may also arise from external disturbances in process operating parameters [11][12][13], in column adsorption capacity, as well as in adsorption isotherm kinetics [14] listed in Table 1.However, for the purpose of this study, only the time-variant variability in u is considered.
In the deterministic robust formulation, the uncertainty set for all disturbances w with L 2 norm less than Γ = σ w u * .Here σ w is the uncertainty level and u * is given by: In order to assess the validity of the linear approximation in Equation ( 13), the robustness of the open-loop system was evaluated by means of Monte Carlo simulations.The result is a probability distribution of the variability of X B (t f ).In each forward simulation, the uncertainty values of the zero-order hold control input, w, were randomly generated from a uniform distribution on [−σ w , σ w ].
The resulting stochastic control signal is thus given by: As outlined in Section 4, two inherently different optimal control strategies are considered and the discrepancy lies in whether Equation ( 11) is augmented with the limit cycle criteria Equation ( 6) or not.Moreover, in order to assess the robustness of the general elution trajectories governed by the optimal control strategy I, these were benchmarked with that of the conventional linear trajectories governed by: where t lin,0 = t 0 + ∆t load + ∆t wash defines the onset of the elution mode and (u lin (t lin,0 ), u lin (t f )) ∈ R 2 are the time-invariant optimization parameters.For the benchmark investigation, an analogous optimal control problem was considered where the quadratic penalty on ∆u in Equation (11a) was removed and bounds on the aforementioned parameters was introduced.In all cases, the first-and second-order spatial derivative of the convection-diffusion Equation ( 1), that governs the dynamics of the mobile phase, have been approximated using a WENO scheme of fifth order with n v = 20 finite volume elements.The number of finite volumes elements is a compromise between accuracy and computational complexity, and is experimentally verified to give adequate representation of the dispersion.Moreover, the temporal horizon was discretized with n e = 2 × 10 2 finite elements with two Radau collocations points in each element and a piecewise constant control discretization with N u = 50 pieces was adopted.Accordingly, the resulting NLP Equation ( 21) has approximately 1.5 × 10 5 NLP variables and was subsequently solved using the primal-dual interior point method IPOPT v.3.11.8 [43] and the linear solver MA57 from HSL [85].The collocation method used corresponds to a fixed-step-size Radau solver.
To verify this temporal discretization, the optimal input and parameters are used to simulate the system using CVode from the SUNDIALS suite [80], which is a variable-step-size, backward-differentiation formula solver with error control.In this study, the robust counter part problem Equation ( 19) was solved on an Ubuntu 12.04 computer with an Intel R Core TM i7-2600 Quad Processor @ 3.40 GHz.The lower-level NLP Equation (21) constrained by the DAE system dynamics was solved using IPOPT to a tolerance of 1 × 10 −13 and a total CPU time of approximately 0.75 h.Moreover, the coupled system of the Lyapunov differential equations and the process model ODEs governed by Equation ( 26) was simulated in the upper-level static optimization problem (Equation (19a)-(19d)) using the sparse representation of the Jacobian together with a connection to the sparse solver in CVode.The total CPU time spent to solve Equation ( 26) was approximately 2.0 h.It is noteworthy that the total CPU time exceeded 48 h (simulating the same system) when sparseness was not considered and when the Jacobian Equation ( 27) was not supplied to the numerical integrator.Hence, this clearly motivates the tools and methods developed for simulating the coupled system of approximately 2.0 × 10 4 state variables (see Figure 5).The remainder of this section is divided into two subsections, where the nominal open-loop and robustified optimal solutions governed by the optimal control strategy I and II are presented in Sections 6.1 and 6.2, respectively.Moreover, the optimization specifications and the associated robustified safety margins for the different values of R, i.e., the penalty on the difference on the zero-order hold control discretization, investigated are summarized in Tables 3 and 4 for comparison purposes.Table 3. Optimization specifications and associated robustified safety margins for the optimal control strategy I and for different values of the penalty R. The optimization specifications correspond to the graphical results depicted in Figures 6 and 10.

Solutions to the Open-loop Optimal
Solutions to the Robust Counterpart Control Problem (Equation ( 11)) Problem (Equation ( 19)) Table 4. Optimization specifications and associated robustified safety margins for the optimal control strategy II and for different values of the penalty R. The optimization specifications correspond to the graphical results depicted in Figures 11 and 14.
Solutions to the Open-loop Optimal Solutions to the Robust Counterpart Control Problem (Equation ( 11)) Problem (Equation ( 19)) The nominal solution (i.e., without the presence of uncertainty) of the open-loop optimal control problem excluding the limit cycle criteria Equation ( 6) is computed over a fixed temporal horizon of ) is the normalized time in column volumes (CV).The specifications in the optimization problem are X B,L = 9.75 × 10 −1 (−), ∆u U = 2.5 × 10 −1 , u L = 5.0 × 10 −2 and u U = 2.0.Moreover, to reduce the state sensitivity of the optimal trajectories to control input uncertainties, we have introduced R as a key specification in the optimization problem.The resulting optimal state and control trajectories are depicted in Figure 6 and the associated optimization specifications are summarized in Table 3.The optimization and verification simulation results are practically identical.The initial control is constrained to be constant over the load and the subsequent wash horizon [t 0 ,t 0 + ∆t load + ∆t wash ] where ∆t wash = 1.0 (CV).Although not considered in the open-loop optimal control problem, the control signal is prescribed a constant value of 1.0 × 10 −1 (moL • kg −1 ) during 2.5 CV regeneration, in order to elute the most retained component, and the initial control, u(t 0 ), is prescribed during the subsequently column re-equilibration period.The optimal state and control trajectories illustrated in Figure 6 show a clear distinction in the system response generated with u and u lin , respectively.As expected from the generic HIC elution mode, u lin and hence the optimal modifier concentration is strictly decreasing in [t lin,0 ,t f ] in order to increase hydrophobicity, see Figure 6d.Governed by the individual component adsorption affinities, components α ∈ {A, B,C} are gradually separated as they traverse the column.It is noteworthy that parameterization Equation (30) for u lin is augmented to the open-loop optimal control problem, however, a zero-order hold fit to that trajectory with N u = 50 segments (without affecting the state trajectories) is considered here for comparison purposes in the robustness assessment.Contrarily, the modifier concentration is freely controlled with the general elution trajectories, u, shown in Figure 6a-c.It is evident that the additional degrees of freedom introduced significantly promotes the recovery yield, Y B (t f ), of the intermediately eluting component B and that Y B (t f ) strictly increases as R → 0. Regardless of R, the solvent strength is gradually increasing until the lower optimal fractionation time, τ 0 , is reached.Moreover, at the onset of the fractionation interval, the slope of the elution trajectory changes sign, and the decreasing solvent strength causes the target component to desorb and ultimately to elute.
It is evident from analyzing Figure 7, where X B is depicted as a function of time for the nominal open-loop controlled trajectories, that the solution of Equation ( 11) yields that the inequality constraint imposed on target component purity is active, i.e., X B (t f ) = X B,L .More importantly, the associated back-off term, augmented to the purity inequality constraint in Equation ( 17), for σ w = 5.0% exhibits that neither of the optimal open-loop controls are robust at the given uncertainty level since Γ C(t f )P (t f )C T (t f ) > 0. This value corresponds to the worst-case purity to be expected and implies that the open-loop controlled system is unable to cope with control input variability experienced during nominal operation.Moreover, the back-off term trajectories pass through a maximum at the onset of the fractionation interval and subsequently decreases as t → t f .As expected, the terminal back-off value decreases with R for the general elution trajectories.It is however noteworthy that the linear elution trajectory is associated with the highest terminal back-off value.This is a key result, as the general elution trajectories not only outperforms the linear in terms of Y B (t f ) but also are inherently more robust.The insight gained from analyzing Figure 7 clearly motivates the scope of this study, and the tools and methods developed here for solving the robust counterpart problem Equation (19).However, since the purity inequality constraint back-off term is determined by linearizing the dynamics around the nominal trajectories, it remains to assess the validity of that approximation.As discussed in Section 4.3, the validity of the robustification strategy is assessed by means of Monte Carlo simulations.A set of stochastic control signals, ũ, given by Equation ( 29) was generated by sampling a uniform distribution on the interval from ±σ w .For each case presented in Figure 6, 5.0 × 10 3 forward simulations were carried out for the uncertainty level of 5.0%.The resulting state and control trajectories are depicted in Figure 8.It is evident that the control signal variability has a direct affect on the individual components' shape and retention time.Most critical is the significant fraction of impurity (component A and C) trajectories spanning inside the fractionation interval and causing degraded target component purity as a result.The corresponding probability distribution of the variability in X B (t f ) is depicted in Figure 9.This figure clearly strengthens the conclusion that the nominal open-loop controlled system is unable to cope with control input variability experienced during nominal operation.More importantly, it is evident that the deterministic back-off value corresponds to a high degree with the worst-case observations of the stochastic distribution of X B (t f ), and consequently, the linearized approximation Equation ( 13) is valid for the uncertainty level considered.Cumulative probability (−) As discussed in Section 1.3, introducing more robustness always give rise to a performance decrease.Hence, the robustification of the open-loop controlled system implies compromising high penalties on the change in the zero-order hold control, yielding optimal solutions with low sensitivity to uncertainty, however, to the expense of an overall lower performance.This is solely the main reason for introducing R as an optimization variable in the upper-level robust counterpart problem Equation (19a)-(19d).It is however noteworthy that neither the stochastic nor the deterministic worst-case purity depicted in Figure 9 is valid for any other open-loop control input associated with a different lower purity requirement than of X B,L = 9.75 × 10 −1 .Therefore, in order to find the optimal robustified controls associated with the nominal open-loop trajectories depicted in Figure 6, an iterative approach was adopted.Hence, this iterative approach implies solving the feasibility problem to the robust counterpart problem Equation (19) when keeping R fixed.The resulting robustified state and control trajectories are depicted in Figure 10 together with their nominal trajectories for comparison purposes.
It is evident from analyzing Figure 10 that the onset of the fractionation interval is slightly shifted towards higher elution volumes for all cases studied.The largest discrepancies between the nominal and the robustified general controls is seen in the temporal region for t > τ 0 .As a consequence, the state trajectory for component A is not affected to a high degree, whereas the state trajectory for the target component B is significantly affected.Moreover, the nominal linear control differ significantly from its robust counterpart both in terms of the initial level and the slope, and consequently, all state trajectories are affected.It is noteworthy that the performance increase in terms of Y B (t f ) for u compared to u lin is even higher for the robustified solutions than for that obtained with the nominal solutions.This can be explained by the lower back-off value associated with the general trajectories.Specifically, the lower purity constraint augmented with the safety margin, X B,L + Γ C(t f )P (t f )C T (t f ), is 9.821 × 10 −1 for R = 0.0 and 9.836 × 10 −1 for u lin .Finally, for the HIC separation case considered in this paper, the target component yield governed by the robustified general elution trajectories is strictly decreasing with the penalty on the change in the zero-order hold control, and consequently, the solution of the robust counterpart problem Equation ( 19) is obtained for R = 0.0, see Figure 10a.This was not expected a priori since the unpenalized input control is associated with the largest purity inequality constraint safety margin.However, the case considered in this study clearly demonstrates the applicability of the tools and methods developed here for formulating and solving the robust counterpart problem.So far, we have shown the optimal open-loop and the robustified elution trajectories for the control strategy excluding the limit cycle criteria.This section is therefore devoted to demonstrate robust optimal control of the comprehensive cyclic-steady-state formalism.As outlined in Section 4.1, computation of CSS solutions over the temporal horizon t ∈ [t 0 ,t f ] requires that open-loop optimal control problem Equation ( 11) is augmented with the limit cycle criteria Equation (6) governing that the state at the initial time is retained at the end of the cycle.The resulting nominal state and control trajectories that conform to the periodic boundary condition over a fixed temporal horizon of ∀ t ∈ [0.0, 25.0] are displayed in Figure 11 and the associated optimization specifications are summarized in Table 4. Here, the initial control, u(t 0 ), is constrained to be constant over the load and the subsequent wash horizon [t 0 ,t 0 +∆t load + ∆t wash ].Moreover, the control in the final re-equilibration horizon [t f − ∆t eq ,t f ], where ∆t eq = 2.0 (CV), is prescribed the constant value of u(t 0 ) through the limit cycle criteria in Equation (6d).In accordance with the optimal open-loop controlled trajectories depicted in Figure 6, the solvent strength is gradually increasing until the lower optimal fractionation time, τ 0 , is reached.For R > 1.0, the slope of the control trajectory changes sign within the fractionation interval in order to prevent the most retained component to desorb.Moreover, at the upper fractionation interval endpoint, the control trajectory drops drastically in order to completely elute component C, and to fulfill Equation (6a), before returning to the initial control level.For these reasons, the upper endpoint of the fractionation interval, τ 0 + ∆τ, do not coincide with the that of the temporal horizon, as was the case for the optimal control trajectories shown in Figure 6.It is also evident from Figure 11 that the penalty on ∆u n significantly affect the state and control trajectories, due to the large difference of the maximum and minimum values of u n , and that the target component recovery yield strictly increases as R → 0.
The deterministic back-off terms determined by linearizing the system dynamics around the nominal trajectories that conform to the CSS criteria are shown in Figure 12.The Γ C(t)P(t)C T (t) trajectories exhibit analogous temporal behavior as for those depicted in Figure 7.Likewise, the terminal back-off value decreases with R. As the purity inequality constraint is active, and terminal back-off value is larger than zero ∀R ∈ [0.0, 1.0], the open-loop controlled system is unable to cope with control input variability experienced during nominal operation.The validity of the deterministic worst-case purity was also here assessed by means of Monte Carlo simulations, as previously outlined in Section 6.1.The probability distribution of the variability in X B (t f ) depicted in Figure 13 is determined by the state trajectories subject to the set of stochastic controls, ũ, shown in Figure 11.The deterministic terminal back-off value corresponds to a high degree with the worst-case observations of the stochastic distribution of X B (t f ) shown in Figure 13, and consequently, the linearized approximation Equation ( 13) is valid for the uncertainty level considered and can therefore safely be exploited in the robustification approach.Cumulative probability (−) The resulting robustified control trajectories, from solving the feasibility problem to the robust counterpart problem Equation (19) when keeping R fixed, are depicted in Figure 14 together with their nominal trajectories for comparison purposes.Again, the largest discrepancy between the nominal and robustified controls are seen in the temporal region around the fractionation interval and τ 0 is slightly shifted towards larger elution volumes, whereas ∆τ is unchanged.It is also evident that the initial control levels almost coincide.Consequently, the state trajectories for component A and C is not affected to a high degree.Hence, in order to meet the robustified purity inequality constraint augmented with the back-off term Equation (17), only the state trajectory of component B is slightly modified.Here, X B,L + Γ C(t f )P (t f )C T (t f ) is [9.815, 9.814, 9.806, 9.803] × 10 −1 for R ∈ [0.0, 0.05, 0.5, 1.0].Hence, there is only a moderate increase in the terminal back-off value as R → 0, and consequently, the solution of the robust counterpart problem Equation ( 19) is obtained for the unpenalized control, see Figure 14a.

Conclusions
External process disturbances are inevitably present and the robustness of the open-loop controlled system with respect to CPP uncertainties, e.g., guaranteeing that QCA requirements are not violated, is of the highest importance in biopharmaceutical process industry.Especially, the solution of the open-loop optimal control problem will yield trajectories that lie on the boundary of the feasible region, and consequently, the open-loop controlled system is unable to cope with control input variability experienced during nominal operation.This paper is therefore concerned with extending the open-loop optimal control problem, constrained by a nonlinear dynamic system governed by PDEs, with its robust counterpart formulation.However, there is a drastic increase in computational complexity while moving from the nominal problem to its robust counterpart.To efficiently and accurately generate optimal elution trajectories in the presence of uncertainty we considered a deterministic robust formulation, that employs a linearization of the uncertainty set, for dynamic systems based on Lyapunov differential equations.While the approximate robust approach is computationally tractable for medium-scale ODE-constrained DOPs, it is not for the original PDE-constrained DOP considered here.One of the main contributions of this paper is therefore the established methodology and the computationally efficient framework for robust optimal control of batch chromatographic separation processes.It is however noteworthy that the generic methods and tools developed here are applicable to any large-scale optimal control problem constrained by PDEs.
Computational tractability was obtained by casting the robust counterpart problem in the frame of bilevel optimal control where the upper level concerns forward simulation of the Lyapunov differential equations and the system dynamics, and the nominal DOP augmented with the robustified inequality constraint margin is considered in the lower level.The lower-level DOP constrained by the system dynamics is transcribed into a sparse NLP using direct collocation, which is then solved by IPOPT supplied with the first-and second-order derivatives of the NLP functions computed using CasADi's algorithmic differentiation while preserving sparsity.The adopted decomposition strategy thereby avoids the full simultaneous discretization of the Lyapunov states, which are too numerous to be treated with direct collocation.However, the coupled system of the Lyapunov differential equations and the process model ODEs significantly increased the problem dimension of the upper level.In order to reduce the over-all execution time it was crucial that the sparsity is used during the simulation and that the structure is used during the assembly of the partial derivatives.This was realized by extending PyFMI to support adding additional equations to the original problem of simulating an FMU as well as to handle the sparse representation of the Jacobian together with the connection to the sparse solver in CVode.Finally, in order to reduce the state sensitivity of the optimal trajectories to control input uncertainties, we have introduced an optimization variable in the upper level that penalizes the change in the zero-order hold control.This parameter influence the smoothness of the control and consequently enhance the robustness, however, to the expense of an overall lower performance.
The advantages of the deterministic robustification methodology were illustrated through the solution of a specific challenging ternary complex mixture separation problem of insulin analogs, with the intermediately eluting component as the target, by hydrophobic interaction chromatography.Moreover, two inherently different optimal control strategies were considered and the discrepancy lies in whether robust counterpart problem is augmented with the cyclic-steady-state criteria or not.The deterministic back-off term, determined by linearizing the system around the nominal open-loop controlled trajectories, for both control strategies showed good agreement with the worst-case observations of the stochastic distribution of purity obtained by means of Monte Carlo simulations.A common trend is that the terminal back-off term value strictly decreased with the value of the penalty on change in the zero-order hold control, hence, yielding solutions with lower sensitivity to control variability.

Figure 2 .
Figure 2. The dynamic optimization framework of the JModelica.orgplatform.The problem is formulated using Modelica and Optimica, which is then transferred to CasADi Interface by the JModelica.orgcompiler.CasADi Interface is a symbolic interface between the compiler, user, and collocation algorithm.After collocation, the problem is finally solved by IPOPT.

Figure 3 .
Figure 3.The simulation capabilities of the JModelica.orgplatform.An FMU is loaded using the package PyFMI, which in turn connects to the simulation package Assimulo which provides access to the solvers.The work-flow is controlled through a user-defined script which additionally allows for providing extra equations to the problem, in this paper the Lyapunov equations.

Figure 4 .
Figure 4.The sparsity pattern of the Jacobian of the spatially discretized, with n v = 20, and ODE transformed system.

Figure 7 .
Figure 7. Target component purity, X B (t), (dashed lines) and associated back-off term, Γ C(t)P(t)C T (t), (solid lines) as a function of time for the nominal elution trajectories depicted in Figure 6.

Figure 8 .
Figure 8. State and control trajectories generated with the set of stochastic control signals, ũ, sampled from a uniform distribution with the uncertainty level σ w = 5.0%.Markers indicate the nominal optimal solution at the Radau collocation points and solid and dashed lines the corresponding simulated response.

Figure 9 .
Figure 9. Probability distribution of the variability in X B (t f ) generated with the set of stochastic control signals, ũ, sampled from a uniform distribution with the uncertainty level σ w = 5.0%.The dashed lines indicate the lower purity constraint, X B,L , and the deterministic worst-case purity, X B,L − Γ C(t f )P (t f )C T (t f ).The solid line indicates the cumulative probability, Φ σ 2 w (X B (t f )).

Figure 10 .
Figure 10.Nominal and robustified state and control trajectories.Markers indicate the robustified solution at the Radau collocation points and solid and dashed lines the corresponding simulated response.The light colored state and control trajectories indicate the nominal solution.

6. 2 .
Optimal Control Strategy II-Nominal and Robustified Elution Trajectories in the Cyclic-Steady-State Operation Mode

Figure 11 .
Figure 11.State and control trajectories generated for X B,L = 9.75 × 10 −1 (−) and R ∈ [0.0, 1.0].The set of stochastic control signals, ũ, was sampled from a uniform distribution with the uncertainty level σ w = 5.0%.Markers indicate nominal state and control trajectories that conform to the cyclic-steady state criteria Equation (6) over the temporal horizon [t 0 ,t f ] at the Radau collocation points and solid and dashed lines the corresponding simulated response.

Figure 12 .
Figure 12.Target component purity, X B (t), (dashed lines) and associated back-off term, Γ C(t)P(t)C T (t), (solid lines) as a function of time for the nominal elution trajectories depicted in Figure11.

Figure 13 .
Figure 13.Probability distribution of the variability in X B (t f ) generated with the set of stochastic control signals, ũ, sampled from a uniform distribution with the uncertainty level σ w = 5.0%.The dashed lines indicate the lower purity constraint, X B,L , and the deterministic worst-case purity, X B,L − Γ C(t f )P (t f )C T (t f ).The solid line indicates the cumulative probability, Φ σ 2 w (X B (t f )).

RRFigure 14 .
Figure 14.Nominal and robust cyclic-steady-state elution profiles for R ∈ [0.0, 1.0].Markers indicate the robustified solution at the Radau collocation points and solid and dashed lines the corresponding simulated response.The light colored state and control trajectories indicate the nominal solution.

Table 1 .
HPLC system component design parameters and HIC column specifics.The adsorption capacity and the self-association parameters listed are equal for all insulin variants α ∈ {A, B,C}.