Flatness-Based Design of Experiments for Model Selection

Proceedings,18th Proceedings,18th IFAC IFAC Symposium Symposium on on System System Identification Identification July 9-11, 9-11, 2018. 2018. Stockholm, Stockholm, Sweden on System Identification Proceedings,18th IFAC Symposium July Sweden Available online at www.sciencedirect.com Proceedings,18th IFAC Symposium July 9-11, 2018. Stockholm, Sweden on System Identification July 9-11, 2018. Stockholm, Sweden

ScienceDirect

IFAC PapersOnLine 51-15 (2018) 233–238

Flatness-Based Design of Experiments Flatness-Based Design of Experiments Flatness-Based Design of Experiments Model Selection Flatness-Based Design of Experiments Model Selection Model Selection Model Selection ∗,∗∗ ∗,∗∗ Moritz Schulze ∗,∗∗ e Schenkendorf ∗,∗∗ ∗,∗∗ Ren´ ∗,∗∗

for for for for

Moritz Schulze ∗,∗∗ Ren´ e Schenkendorf ∗,∗∗ Moritz Schulze ∗,∗∗ Ren´ e Schenkendorf ∗,∗∗ Moritz Schulze Ren´ e Systems Schenkendorf ∗ ∗ of Energy and Process Engineering, TU ∗ Institute ∗ Institute of Energy and Process Systems Engineering, TU Institute of Energy and Process Systems Engineering, TU Braunschweig, Franz-Liszt-Straße 35, 38106 Braunschweig (Germany), Braunschweig, Franz-Liszt-Straße 35, 38106 Braunschweig (Germany), ∗ Institute of Energy and Process Systems Engineering, TU Braunschweig, Franz-Liszt-Straße 35, 38106 Braunschweig (Germany), (e-mail: [email protected]). (e-mail: [email protected]). ∗∗ Braunschweig, Franz-Liszt-Straße 35, 38106(PVZ), Braunschweig (Germany), (e-mail: [email protected]). ∗∗ of Pharmaceutical Engineering TU Braunschweig, ∗∗ Center of(e-mail: Pharmaceutical Engineering (PVZ), TU Braunschweig, ∗∗ Center [email protected]). Center of Pharmaceutical Engineering (PVZ), TU Braunschweig, Franz-Liszt-Straße 35a, 38106 (Germany). Franz-Liszt-Straße 35a, 38106 Braunschweig Braunschweig (Germany). ∗∗ Center of Pharmaceutical (PVZ), TU Braunschweig, Franz-Liszt-Straße 35a, Engineering 38106 Braunschweig (Germany). Franz-Liszt-Straße 35a, 38106 Braunschweig (Germany). Abstract: High High product product standards standards and and cost-efficient cost-efficient robust robust process process designs designs are are key key issues issues in in Abstract: Abstract: High product standards and cost-efficient robust understanding process designsof are key issues in the pharmaceutical industry and require thorough system the underlying the pharmaceutical industry and require thorough system understanding of the underlying Abstract: High product standards and cost-efficient robust process designsofto are key issues in the pharmaceutical industry and require thorough understanding the underlying physical phenomena. Model-based approaches have system been proven favorable gain valuable physical phenomena. Model-based approaches have been proven favorable tothe gain valuable the pharmaceutical industry and require thorough system understanding of underlying physical phenomena. Model-based approaches been proven to gain valuable system insights. insights. However, qualitative inferences have can be drawn only favorable if the the model model is adequately system However, qualitative inferences can drawn only if is physical phenomena. Model-based have been proven to gain system However, inferencesdesigned can be be drawn only favorable ifIn thethis model is adequately adequately derived insights. and validated by – qualitative at best best – –approaches optimally designed experiments. In this paper, wevaluable present derived and validated by – at optimally experiments. paper, we present system insights. However, inferences can design be drawn only ifIn thethis model is adequately derived and validated by – qualitative at best – optimally designed experiments. paper, weproblems present a new model inversion strategy to solve model-based of experiments (MBDoE) a new model inversion strategy to solve model-based design of experiments (MBDoE) problems derived and validated – maximize at best – optimally designed experiments. this paper,byweproblems present afornew model inversion strategy to solve of experiments (MBDoE) model selection; i.e.,by we themodel-based difference ofdesign competing modelIncandidates candidates using the for model selection; i.e., we maximize the difference of competing model by problems using the a new model inversion strategy to solve model-based design of experiments (MBDoE) for model selection; i.e., we maximize the difference of competing modeltheir candidates by using the differential flatness concept. Differentially flat systems can represent states and controls differential flatness i.e., concept. Differentially flat systems can represent states and controls for model selection; we maximize theand difference of competing modeltheir candidates by using differential flatness concept. Differentially flat can By represent their states and controls analytically by the the so-called so-called flat output its systems derivatives. By maximizing the difference of the the analytically by flat output and its derivatives. maximizing the difference of the differential flatness concept. flat Differentially flat systems can By represent their the states and controls analytically by the so-called output and its derivatives. maximizing difference of the flat outputs of all considered model candidates, we efficiently obtain analytical model-revealing flat outputs of all considered model candidates, we efficiently obtain analytical model-revealing analytically by theconsidered so-called flat output and itsthese derivatives. Byobtain maximizing theare difference of the flat outputs of all model candidates, we efficiently analytical model-revealing control actions for the MBDoE problem. When optimized control actions implemented, control actions for the MBDoE problem. When these optimized control actions are implemented, flat outputs of all considered model candidates, we efficiently obtain analytical model-revealing control actionsexperimental for the MBDoE When to these optimized controlmodels actionsmore are implemented, the resulting resulting dataproblem. are expected expected invalidate incorrect reliably. We We the experimental data are to invalidate incorrect models reliably. control actions foreffectiveness the MBDoE WhenMBDoE these optimized control actionsmore are implemented, the resulting experimental data are proposed expected to invalidate incorrect models more reliably. We demonstrate the ofproblem. the proposed concept with potential potential model candidates demonstrate the effectiveness of the MBDoE concept with model candidates the resulting the experimental data areofproposed expected to invalidate incorrect models more We demonstrate effectiveness of the MBDoE concept with potential modelreliably. candidates that realizations a reaction system. that describe describe different different realizations ofproposed a mass mass action action reaction system. demonstrate effectiveness of theof MBDoE concept with potential model candidates that describethe different realizations a mass action reaction system. © 2018, IFAC (International Federationofof aAutomatic Control) Hosting by Elsevier Ltd. All rights reserved. that describe different realizations mass action reaction system. Keywords: Keywords: model model selection, selection, differential differential flatness, flatness, design design of of experiments, experiments, optimal optimal control, control, Akaike Akaike Keywords: model selection, differential flatness, design of experiments, optimal control, Akaike information criterion information criterion Keywords: selection, differential flatness, design of experiments, optimal control, Akaike informationmodel criterion information criterion 1. INTRODUCTION called flat flat system outputs outputs and their their derivatives. derivatives. Thus, Thus, by by 1. called 1. INTRODUCTION INTRODUCTION called flat system system outputs and and their strategy derivatives. Thus, by implementing a differential flatness for MBDoE, implementing a differential flatness strategy for MBDoE, 1. INTRODUCTION called flat system outputs and their strategy derivatives. by implementing a differential flatness for Thus, MBDoE, we can avoid simplifying assumptions and the need to we can avoid simplifying assumptions and the need to solve solve Profit margins in the highly competitive and regulated implementing a differential flatness strategy for MBDoE, can avoidequations simplifying assumptions and the need to solve differential numerically. Profit margins in the highly competitive and regulated we differential equations numerically. Profit margins in the highly and regulated pharmaceutical industry have competitive been continually continually declining we can avoidequations simplifying assumptions and the need to solve differential numerically. pharmaceutical industry have been declining Profit margins in the highly competitive and regulated pharmaceutical industry have been continually declining The theory of differential flatness over the last two decades because of increasing research differential equations numerically. The theory of differential flatness for for dynamic dynamic systems systems was was over the last two decadeshave because ofcontinually increasing declining research The theory of differential flatness pharmaceutical industry been for dynamic systems was derived by Fliess et al. (1992). Since then, flatness-based over the last two decades because of increasing research and development costs (Behr et al., 2004). For instance, and development costs (Behr et al., 2004). For instance, derived by Fliess et al. (1992). Since then, flatness-based theory of differential flatness forresearch dynamic systems was over thethe last two process decades because of2004). increasing research derived by Fliess et widely al. (1992). then, flatness-based methods have been usedSince in and industry and development costs (Behr etofal., For instance, during whole chain producing active phar- The methods have been widely used in research and industry during the whole process chain ofal., producing active pharderived by Fliess et al. (1992). Since then, flatness-based and development costs (Behr et 2004). For instance, methods have been widely used in research and industry primarily in in controller controller design design and and trajectory trajectory planning planning of of during the whole process(APIs), chain of producing activeunderphar- –– primarily maceutical ingredients thorough system maceutical ingredients thorough system undermethods have beensystems widely usedand in and research and industry during theofwhole process(APIs), chain of producing active phar–electromechanical primarily in controller design trajectory planning of maceutical ingredients (APIs), thorough system under(Franke Robenack, 2013). standing the underlying physical processes is key to electromechanical systems (Franke and Robenack, 2013). standing of ingredients the underlying physical processes is key to primarily inMBDoE, controller design and trajectory planning of maceutical (APIs), thorough system undersystems (Franke and 2013). Considering strategies based onRobenack, differential flatstanding the product underlying physical processes is key to –electromechanical meet high standards and to costmeet the the of high product standards and to enable enable costConsidering MBDoE, strategies based on differential flatelectromechanical systems (Franke and Robenack, 2013). standing of the underlying physical processes is key to MBDoE, based on differential ness have have rarely rarely beenstrategies considered (Schenkendorf et flatal., meet thedrug highdevelopment product standards and to enable efficient and manufacturing. manufacturing. To costgain Considering ness been considered (Schenkendorf et al., efficient drug development and To gain MBDoE, strategies based on differential meet thedrug high productthestandards andoftofirst-principles enable ness have rarely been considered (Schenkendorf et flatal., 2012; Schenkendorf and Mangold, 2014) and to our knowlefficient development and manufacturing. To costgain Considering true system insights, application 2012; have Schenkendorf and Mangold, 2014) and to our knowltrue system insights, the application of rarely been considered (Schenkendorf et al., efficient drugbeen development and manufacturing. To gain 2012; Mangold, 2014) and to our knowltrue system insights, thebeneficial. application of first-principles first-principles edge have applied to selection problems models has proven beneficial. Process systems en- ness edge Schenkendorf have not not been beenand applied to model model selection problems models has been proven Process systems en2012; Schenkendorf and Mangold, 2014) and to our knowltrue system insights, the application ofof first-principles edge have not been applied to model selection problems in the literature. models has been proven beneficial. Process systems engineers, however, often face a plurality mathematical gineers, has however, often face a plurality of mathematical in the literature. not been applied to model selection problems models been proven beneficial. systems in the have literature. gineers, however, often a plurality of mathematical model candidates candidates which face describe theProcess reaction networkenof edge model which describe the reaction network of In this paper, we in the literature. gineers, however, often face a plurality of mathematical In this paper, we extend extend the the application application of of the the flatness flatness model candidates the reaction network of the process process under which study describe equally well. well. In addition addition to the the In this paper, we extend the application of flatness the under study equally In to approach within the MBDoE framework to the class of model candidates which describe the reaction network of approach within the MBDoE framework to the class of the process under study equally well. In addition to the parameter estimation estimation problem, problem, the the identification identification of of the In thisaction paper, wethe extend the application of the the class flatness parameter of approach within MBDoE framework to mass systems because of their fundamental imthe processestimation under study equallyisthe well.identification In for addition of to the the mass action systems because of their fundamental imparameter problem, most model approach within the MBDoE framework to the classimof most suitable suitable model candidate candidate is crucial crucial for model-based model-based action systems because networks of their in fundamental pact on on modeling of reaction reaction networks in (bio)chemistry parameter estimation problem, identification of the mass pact modeling of (bio)chemistry most suitable model By candidate isthe crucial for model-based analysis and design. following the model-based design mass action systems because of their fundamental impact on modeling of reaction networks in (bio)chemistry analysis and design. By following the model-based design (Voit et al., 2015). We are aiming at falsifying improper most suitable model candidate is crucial fora model-based (Voit on et modeling al., 2015).ofWe are aiming at falsifying improper analysis and design. By following thesolving model-based of experiments experiments (MBDoE) principle, modeldesign iden- pact reaction networks in (bio)chemistry (Voit al., 2015). are aiming at falsifying improper of (MBDoE) principle, solving a model idenmodelsetfrom from a set set of ofWe competing model candidates in two two analysis and design. By following the model-based design models a competing model candidates in of experiments (MBDoE) principle, solving a model identification problem problem for for dynamic dynamic systems systems requires requires optimal optimal (Voit et al., 2015). We are aiming at falsifying improper tification models from a set of competing model candidates in two steps. First, we show that the individual models satisfy of experiments (MBDoE) principle, solving a model iden- steps. First, we show that the individual models satisfy tification problem for dynamic requires control Commonly used collocation and control fromflatness a set ofconditions, competing model candidates in two control actions. actions. Commonly usedsystems collocation and optimal control models steps. First, we show that the individual models satisfy differential and second, we solve solve an tification problem for dynamic systems requires optimal differential flatness conditions, and second, we an control actions. Commonly used collocation and control steps. parameterization techniques result in complex optimizaFirst,flatness we show that uses the individual models satisfy differential conditions, and second, we solve an parameterization techniques result in complex optimizaoptimization problem that the flatness property. We control actions. Commonly used collocation and control optimization problem that uses the flatness property. We parameterization techniques result in complex optimization problems problems that that require require substantial substantial computational computational time time differential flatness conditions, and second, we solve We an optimization problem that uses the flatness tion efficiently derive derive analytical optimal controlproperty. actions for parameterization result in computational complex optimizaefficiently analytical optimal control actions for tion problemssolvers thattechniques require substantial time and efficient (Biral et al., 2016). Alternatively, the optimization problem that uses the flatness property. We and efficient solvers (Biral et al., 2016). Alternatively, the efficiently derive analytical optimal control actions for subsequent validation validation experiments experiments to to figure figure out out the the most most tion problems that require computational time subsequent and efficient solvers (Biral substantial etenables al., 2016). Alternatively, the efficiently differential flatness concept effective model selecselecderive analytical optimal control actions for subsequent validation to figure out the most differential flatness concept enables effective model suitable model from aaexperiments set of potential model candidates and efficient solvers (Biral et al., 2016). Alternatively, the suitable model from set of potential model candidates differential conceptframework. enables effective modeltheory, selec- subsequent validation experiments to figure out the most tion within within flatness the MBDoE MBDoE In systems systems suitable model from a set of potential model candidates tion the framework. In theory, describing different realizations of a mass action reaction differential concept enables effective model selection within flatness the MBDoE framework. In systems theory, describing different realizations of a mass action reaction differential flatness ensures that the the system states and suitable fromrealizations a set of potential model candidates describing different of a mass action reaction differential flatness ensures that system states and system. model tion within functions the MBDoE framework. In systems theory, system. differential flatness ensures that the system states and the control can be expressed analytically by sodifferent realizations of a mass action reaction the control functions can be expressed analytically byand so- describing system. differential flatness ensures that the system statesby the control functions can be expressed analytically so- system. the control functions can be expressed analytically so- Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2018, IFAC (International Federation of AutomaticbyControl)

Copyright © 2018 IFAC 233 Copyright 2018 responsibility IFAC 233Control. Peer review©under of International Federation of Automatic Copyright © 2018 IFAC 233 10.1016/j.ifacol.2018.09.140 Copyright © 2018 IFAC 233

2018 IFAC SYSID 234 July 9-11, 2018. Stockholm, Sweden

Moritz Schulze et al. / IFAC PapersOnLine 51-15 (2018) 233–238

2. METHODS Dynamic system behavior is often mathematically described by ordinary differential equations (ODEs) which read in their state-space representation as dx = f (x(t), u(t), θ), (1) dt where the dynamic of the system is a function of the dynamic states x(t) ∈ Rn , the system input u(t) ∈ Rm , and the parameter vector θ ∈ Rp . The non-linear function f : R n × Rm × R p − → Rn represents the vector field of the dynamic system. The corresponding output equation is y(t) = h(x(t), u(t), θ) (2) q n m p q → R as the with outputs y(t) ∈ R and h: R × R × R − (non-)linear output function. For readability reasons, the notion of time dependency and the parameter vector are omitted but are specified when the context demands it.

flat outputs (3) followed by a subsequent validation step in the form of equations (4), (5), and (6) has proven favorable. This trial-and-error approach is supported by the fact that many technical systems are flat systems, and that often, informative flat outputs have physical meaning, similar to Lyapunov functions (Adamy, 2014). 2.2 Design of experiments for model selection In the literature one can distinguish between statistical and model-based design of experiments, where methods of the former, e.g., response surface methods, are popular and widely used because of their simplicity but are not able to properly cope with complex dynamic systems for the same reason (Franceschini and Macchietto, 2008). A generic scheme of the proposed MBDoE strategy for model selection is given in Fig. 1. Typically, MBDoE results in a dynamic optimization problem of the form max J (φ) φ∈Φ

2.1 Differentially flat systems

(7) geq (φ) = 0, gineq (φ) ≤ 0. Here, J is a scalar objective function, and φ(u, x0 ) is the design vector within the design space Φ specified by the control u and the initial states x0 = x(t = 0). The optimization problem is subject to equality constraints geq and inequality constraints gineq which, for instance, might be introduced to limit the concentration of a species or to ensure positive values of the system states. Various scalar objective functions for MBDoE have been proposed in the literature; an overview can be found in Verheijen (2003). Thus, an optimal control problem is given, where the continuous time-dependent control inputs u(t) are parameterized by solving the underlying ODE system (1) numerically while aiming to maximize some distance measure between the competing model candidates. s.t.

A non-linear dynamic system {(1),(2)} is called differentially flat or shortly flat if there exists an output vector ˙ . . . , u(s) ), (3) yflat = hflat (x, u, u, with a finite value s ∈ N, referred to as a flat output which fulfills the following conditions: (1) The states and the controls can be described as a function of the flat output and its derivatives: x = Ψx (yflat , y˙ flat , . . . , yflat

(r)

),

(4)

flat (r+1)

). (5) u = Ψu (yflat , y˙ flat , . . . , y (2) The dimensions of the control and the flat output vector are equal: dim yflat = dim u. (6) In general, r ∈ N is not known a priori, except for singleinput and single-output (SISO) systems, where r = n − 1 (Adamy, 2014). There is an infinite number of flat outputs, and often, they are a function of the states only (Adamy, 2014), where valuable information about the flat output design can be extracted from graph theory (Wey, 2002; Richard et al., 2002). Flat outputs might be fictitious or real, i.e., the flat outputs have no physical meaning, or they might be equal to measurable quantities. In the latter case, the controls are also referred to as flat inputs (Waldherr and Zeitz, 2008). Flatness-based methods are also being researched and used for partial differential equations (PDEs); see Wagner and Meurer (2018) and references therein. In practice, there is no general method for either determining if a system is flat or constructing a flat output. Several exceptions exist: Systems that are linearizable by static feedback control are flat systems by definition. In the case of linear systems, differential flatness is equal to the system being controllable and vice versa (Rigatos, 2015). The construction of flat inputs was shown for SISO systems and a limited set of multi-input and multi-output (MIMO) cases (Waldherr and Zeitz, 2010). Further methods of determining flat outputs are currently being researched (Kolar et al., 2016; Franke and Robenack, 2013; Victor et al., 2014). In practice, an expert guessing strategy to identify 234

The proposed flatness-based approach for MBDoE performs differently. It operates in an inverse manner. The difference between the model outputs is maximized by optimized flat output curves; i.e., the corresponding controls u and the initial states x0 are recalculated and have to be equal for all considered model candidates and their optimized flat outputs. Assuming M model candidates and taking the Euclidean metric as distance measure, the optimization problem (7) reads as tf M −1  M   (r) max yi (yflat,i , y˙ flat,i , . . . , yflat,i ) yflat

t0

i=1 j=i+1

− yj (yflat,j , y˙ flat,j , . . . , yflat,j s.t.

∆u(yflat ) = 0,

(r)

)

2

dt (8)

∆x0 (yflat ) = 0, geq (yflat ) = 0, gineq (yflat ) ≤ 0. The inputs and outputs are expressed according to equations (2), (4), and (5). The Delta function ∆(·) measures the differences between the recalculated inputs and the initial states for all M model candidates. The flat outputs themselves are specified by time-dependent empirical

2018 IFAC SYSID July 9-11, 2018. Stockholm, Sweden

Moritz Schulze et al. / IFAC PapersOnLine 51-15 (2018) 233–238

235

start >1 Model building

DoE

Experiment

Parameter estimation

Number of models left?

Model rejection

1 end

0

Fig. 1. MBDoE work-flow for model selection basis functions yflat,i = yflat,i (θ,i , t) ∀ i ∈ {1, ..., M } and their meta-parameters θ,i . As the inputs and the outputs are available in closed form, the optimal control problem given in (7) is transformed in an algebraic nonlinear optimization problem. The gained optimal control profiles and the initial conditions are subsequently used to run a new experiment which is expected to provide more informative data in terms of model falsification. After the experiment is run and the parameters estimated, model invalidation is performed. Here, we can distinguish between several approaches for model selection: testing strategies, goodness-of-fit measures, cross-validation, and criteria based on information theory (Burnham et al., 2002). A well-accepted metric based on information theory is the corrected Akaike information criterion (AIC) 2K(K + 1) AICc = −2 log L(θ) + 2K + , (9) N −K −1 where K is the total number of estimated regression parameters, N is the number of data sample points, and L is the likelihood function, which measures the mismatch between simulation results and the measurement data. Please note, that the use of the corrected AIC is recommended for values N/K below 40 (Burnham et al., 2002). 2.3 Basis functions

D, are produced. To describe this process mathematically, three different potential model candidates are proposed. The candidates were deliberately chosen to integrate typical phenomena that are considered in model selection when dealing with chemical reaction systems like adding/removing of intermediates and inclusion of back reactions. Product E is formed by: 1) reactant B and intermediate D, 2) intermediate D solely, and 3) intermediate D but the reaction is reversible. The described reaction mechanisms are summarized and shown below. k

5 A + 2 B −−1→ C k

C −−2→ A + D

B + D −−→ E k

5 A + 2 B −−1→ C k

C −−2→ A + D

D −−→ E k

5 A + 2 B −−1→ C k

3. CASE STUDY The studied mass action system forms product E from two reactants A and B while two intermediates, C and 235

Model 2

k3

C −−2→ A + D

Selecting suitable, preferably analytic basis functions is a non-trivial task because of the complex expressions of the inverse system’s states and controls resulting from the nonlinear structure of the underlying system of differential equations. Often, higher-order derivatives appear, where small changes in these derivatives may cause large effects in the states and controls. Additionally, differentiability and the end conditions have to be taken care of, notably the initial conditions of the states, which are not part of the flat output vector. Previously free to choose when the system of differential equations was numerically integrated, the initial conditions are now bounded by the flat output and its derivatives. Further constraints, e.g., lower and upper limits on states and constraints or restrictions on operational characteristics, can be incorporated as well, but have to be considered in parallel when choosing the basis functions. Concerning model selection, where constraints are integrated into the optimization problem, an evolutionary approach starting from simple functions with few meta-parameters as lower order polynomials towards functions, that ensure considerable design freedom like splines, is recommended.

Model 1

k3

Model 3

k3

−− D − −E k4

Model 1 is assumed to represent the true model and is used to generate simulated data, i.e., in-silico measurement data with a zero intercept and a constant variance. In principle, however, the reactants and products could be measured in the laboratory, and both reactants might be added independently via pumps and used as control inputs for MBDoE. Thus, for example, the state equations of the dynamic system of Model 1 according to equation (1) are 2/5

x˙ 1 = −k1 x1 x2 + k2 x3 + u1 2 2/5 x˙ 2 = − k1 x1 x2 − k3 x2 x4 + u2 5 1 2/5 x˙ 3 = k1 x1 x2 − k2 x3 5 x˙ 4 = k2 x3 − k3 x2 x4 x˙ 5 = k3 x2 x4 .

(10)

The state vector x comprises the concentrations of the five participating species A to E. Typically, not all quantities T are measurable. Thus, we assume y1 = (x1 x5 ) as the output vector for Model 1. Please note that measuring x1 does not result in a flat output but might be reconstructed by an alternative flat output as shown below. Model 2 and Model 3 can be assembled analogically.

2018 IFAC SYSID 236 9-11, 2018. Stockholm, Sweden July

Moritz Schulze et al. / IFAC PapersOnLine 51-15 (2018) 233–238

4.1 Flat output vector As x5 does not appear on the right side of the dynamic equations in (10), it must be included in the flat output vector. A self-evident guess for additional flat outputs with a physical meaning are the concentrations themselves. For instance, a flat output candidate for Model 1 is yflat,1 = T (x4 x5 ) ; i.e., the conditions (4), (5), and (6) are fulfilled. The resulting algebraic terms of the states and controls are not shown here for the lack of space. Two situations occur where the system is not controllable with the specified flat output, i.e., when x4 = 0 or x˙ 5 = 0. By changing the flat output to   x2 flat,1 y = , (11) x5 we circumvent the singularity at x˙ 5 = 0; this effect is known in the literature as an extrinsic singularity (Kaminski et al., 2018). In either case, the singularities lie on constraint borders of the optimization problem if assumed that the concentrations can be interpreted only when they are positive. Accordingly, it can be shown that the output vector (11) is a flat output for Model 2 and Model 3 alike. Thus, all of the model candidates’ dynamic systems are differentially flat. 4.2 Trajectory planning As mentioned previously, the flatness-based MBDoE strategy benefits from the easiness in system trajectory planning, i.e., for desired system dynamics the related control actions are immediately available. This is in particular true, when the flat output is physically linked to the system. Therefore, in the case study at hand, a target trajectory y∗ (t) = yflat,∗ (t) can be designed by choosing basis functions that express the desired time-varying concentration profiles. In a first step before model selection, we aim at transferring the system to a stable operational point, where, for example, Hagenmeyer and Zeitz (2004) used polynomial functions of degree five to realize this transition. In contrast, in the proposed reaction system the polynomial approach cannot be used due to the mentioned sensitivity effect induced by the used derivatives. In higher-order polynomials, oscillating effects, also known as Runge’s phenomenon, occur and may cause difficultto-implement control actions. In all three models, fourthorder derivatives of the flat output appear in the controls. Simulations have shown that the states and controls are very sensitive to minor changes in these derivatives. At the same time, both model outputs of the flat output vector have to be analyzed in parallel; i.e., the two presumably autonomous design functions cannot be chosen independently if the resulting behavior of the system is restricted to bounds. In the following, we consider Model 1 and illustrate how to design the desired system trajectory and how it relates to control actions. The concentration of reactant B takes the role of a constant substrate and is therefore expressed as x∗2 (t) = const. (12) 236

Controls u

4. RESULTS

1.5

u1 u2

1 0.5 0 0

100

200 300 Time t

400

500

Fig. 2. Control actions needed for the target trajectories The transient concentration profile of product E is assumed to gradually reach a plateau of desired product concentration, at which, for instance, the product solution could be removed from the reactor. We can realize the desired behavior using exponential functions equal to    T +t t g(t) = R 1 − exp − , (13) T T with constants R and T . From the analytic expressions of the intermediate concentrations x3 (t) and x4 (t) of the inverse system, we obtain that the first and second derivatives of the concentration of product E x5 (t) at the initial point t0 = 0 have to vanish. Otherwise, intermediates C and D have to be available in stock in order to add them to the reactor at the beginning of the reaction. Therefore, we expand the formula in (13) by adding another summand and are left with    t T1 + t ∗ x5 (t) = R1 1 − exp − T1 T   1  (14) T2 + t t , + R2 1 − exp − T2 T2 and the side condition R2 R1 + 2 = 0. (15) T12 T2 If a substrate concentration of 15 and a product concentration of around 20 is targeted, the constants in equation (14) can be chosen according to the values shown in Table 1. The remaining constant is calculated from the side condition in (15). Results for the controls and states Table 1. Constants for the target trajectory Constant R1 T1 T2

Value 35 60 40

are shown in Fig. 2 and Fig. 3, respectively. The lines represent the solutions obtained from the inverse Model 1, i.e., by applying the differential flatness concept. The solutions of the ODE system, in turn, were numerically integrated using the controls obtained from the inverse system, and are shown as dots in Fig. 3. As expected, the values coincide perfectly. 4.3 Model selection A translation of the non-linear programming problem (8) into an implementation-related optimization problem results in

Moritz Schulze et al. / IFAC PapersOnLine 51-15 (2018) 233–238

States x

30

x1 x3 x5

x2 x4 ODE

Controls u

2018 IFAC SYSID July 9-11, 2018. Stockholm, Sweden

20

0.8 0.6 0.4 0.2 0

237

u1 u2

0

100

10

200 300 Time t

400

500

Fig. 5. Optimal experimental conditions 0

0

100

200 300 Time t

400

500

Fig. 3. States of the target trajectory max θ

n −1 t −1 M 

M  

n −1 t −1 M 

M 

k=0 i=1 j=i+1

s.t.

k=0 i=1 j=i+1

yi (θi , tk ) − yj (θj , tk )

2

(ui (θi , tk ) − uj (θj , tk ))2 = 0

(θj+1 , t0 ), i ∈ {1, 2}, xji (θj , t0 ) = xj+1 i j ∈ {1, . . . , M − 1}, xi (θ, t0 ) = 0, i ∈ {3, 4, 5}, xi (θ, tk ) ≥ 0, i ∈ {1, . . . , 5}, k ∈ {0, . . . , nt − 1}, xi (θ, tk ) ≤ 50, i ∈ {1, . . . , 5}, k ∈ {0, . . . , nt − 1}, ui (θ, tk ) ≥ 0, i ∈ {1, 2}, k ∈ {0, . . . , nt − 1}, (16) with nt = 30 the number of data time points. According to the inequality constraints on the concentration states, lower and upper bounds are additionally introduced for the design variables. It was assumed that one experiment was run before the optimization. The model parameters were regressed with the results for the corrected AIC shown in the second column in Table 2. The corresponding results of Table 2. AICc Model 1 2 3

not optimized −170.14 −110.18 −141.20

optimized −350.76 +193.31 +186.51

the model outputs are shown in Fig. 4. As can be observed in the plot, all of the models are able to accurately predict their outputs for the first data set. Model selection is not possible at this stage.

Outputs y

Benchmark case. The basis functions for the flat output vector (11) were taken as first-order polynomials resulting y1 (x1 ) y1 (x5 )

40

y2 (x1 ) y2 (x5 )

y3 (x1 ) y3 (x5 )

20 0

0

100

200 300 Time t

400

500

in a 12-dimensional design space. If the initial equality constraints on species C, D, and E are neglected, the solver converges to the objective function value 1.5e5. It circumvents the dynamics of the underlying systems ending up in a situation where only the flat output is nonzero. Thus, the objective function is maximized when the model output of one model is zero, while the model outputs of the two remaining models are at the objective function’s maximum respecting the bound constraints given in (16). In order to ensure an acceptable signal-to-noise ratio, another inequality constraint was introduced: n t −1 xi (tk ) ≥ nt α, i ∈ 1, . . . , 5. (17) k=0

Here, α is a design constant and represents an artificial limit, that depends on the measurement noise of the used measuring device and was set to 5 in the case at hand.

Global optimization using splines. To circumvent difficulties when choosing basis functions, as described in the trajectory planning section, and to omit limitations on the profile of the concentration states, we turn to piecewise functions, because they allow significant freedom when considering the shape of the curves. Commonly used piecewise functions are polynomials within the context of Bsplines. Extensive information about B-splines and splines in general is given in the book by de Boor (2001). In the upcoming optimization problem, all of the flat outputs’ basis functions were expressed as B-spline curves. The number of control points of the B-spline curves, i.e., the number of meta-parameters, was chosen to be 18, leading to a 54-dimensional design space. The subsequent global optimization was done using a multi-start approach in combination with a gradient-based solver. The resulting optimal control is shown in Fig. 5. Following optimization, new in-silico measurement data were created. Together with the first data set, the outputs of the models were regressed anew. The outcome is illustrated in Fig. 6 and the corrected AIC is given in the third column in Table 2. It clearly can be observed from the plot and the AICc values that both data sets can be accurately reproduced solely by Model 1. Additionally, Model 3 performs slightly better than but qualitatively equal to Model 2, as Model 3 is a super-structure model of Model 2. In conclusion, Model 2 and Model 3 can be rejected after the run of a single additional experiment, and Model 1 is selected as the most promising candidate. 5. CONCLUSION We have illustrated a novel approach to the model-based design of experiments for model selection using differen-

Fig. 4. Results for the three models before optimization 237

2018 IFAC SYSID 238 9-11, 2018. Stockholm, Sweden July

Moritz Schulze et al. / IFAC PapersOnLine 51-15 (2018) 233–238

30 data

Outputs y

y (x1 ) ydata (x5 )

1

y (x1 ) y1 (x5 )

2

y (x1 ) y2 (x5 )

3

y (x1 ) y3 (x5 )

20

10

0

0

100

200 300 Time t

400

500

Fig. 6. Model outputs after MBDoE tially flat methods. Compared to commonly used methods, our approach circumvents the numerical approximations needed for solving the emerging optimal control problem, because in differentially flat systems, the states and controls can be derived analytically without integrating the underlying system of differential equations. Moreover, a supply of analytic gradients to the optimization solver, even for highly non-linear systems, is possible. We successfully applied the strategy to a non-linear mass action system case using B-splines as the parameterization technique, thus obtaining optimal experimental conditions for the subsequent experiment. We successfully falsified improper models and selected the most promising one after a single additional experiment. Furthermore, we showed how the differentially flat methods can be used at the same time to design feedforward control for the desired system trajectory of the study case. Moreover, we also stressed that the determination if a non-linear system is flat and the computation of flat outputs are non-trivial tasks and subject of ongoing research. In future, we plan to apply the proposed flatness-based design of experiments strategy to a reaction system of pharmaceutical relevance and to validate the results with laboratory experiments. Furthermore, we will also focus on systematic ways of finding flat outputs and the subsequent choice of basis function while including noisy observations and uncertainties of initial conditions and model parameters, too. REFERENCES Adamy, J. (2014). Nichtlineare Systeme und Regelungen. Springer Vieweg, Berlin, Heidelberg. Behr, A., Brehme, V., Ewers, C., Gr¨ on, H., Kimmel, T., K¨ uppers, S., and Symietz, I. (2004). New Developments in Chemical Engineering for the Production of Drug Substances. Engineering in Life Sciences, 4(1), 15–24. Biral, F., Bertolazzi, E., and Bosetti, P. (2016). Notes on Numerical Methods for Solving Optimal Control Problems. IEEJ Journal of Industry Applications, 5(2), 154–166. Burnham, K.P., Anderson, D.R., and Burnham, K.P. (2002). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Springer, New York, 2nd edition. de Boor, C. (2001). A Practical Guide to Splines. Number 27 in Applied Mathematical Sciences. SpringerVerlag, New York. 238

Fliess, M., L´evine, J., Martin, P., and Rouchon, P. (1992). Sur les syst`emes non lin´eaires diff´erentiellement plats. CR Acad. Sci. Paris, 619–624. Franceschini, G. and Macchietto, S. (2008). Model-Based Design of Experiments for Parameter Precision: State of the Art. Chemical Engineering Science, 63(19), 4846– 4872. Franke, M. and Robenack, K. (2013). On the Computation of Flat Outputs for Nonlinear Control Systems. In European Control Conference (ECC), 167–172. IEEE. Hagenmeyer, V. and Zeitz, M. (2004). Flachheitsbasierter Entwurf von Linearen Und Nichtlinearen Vorsteuerungen (Flatness-Based Design of Linear and Nonlinear Feedforward Controls). at – Automatisierungstechnik/Methoden und Anwendungen der Steuerungs-, Regelungs- und Informationstechnik, 52(1), 3–12. Kaminski, Y.J., L´evine, J., and Ollivier, F. (2018). Intrinsic and apparent singularities in differentially flat systems, and application to global motion planning. Systems & Control Letters, 113, 117–124. ¨ and Kolar, B., Kaldm¨ae, A., Sch¨oberl, M., Kotta, U., Schlacher, K. (2016). Construction of Flat Outputs of Nonlinear Discrete-Time Systems in a Geometric and an Algebraic Framework. IFAC-PapersOnLine, 49(18), 796–801. Richard, P., Buisson, J., and Cormerais, H. (2002). Analysis of Flatness using Bond Graphs and Bicausality. IFAC Proceedings Volumes, 35(1), 25–30. Rigatos, G.G. (2015). Nonlinear Control and Filtering Using Differential Flatness Approaches, volume 25 of Studies in Systems, Decision and Control. Springer International Publishing, Cham. Schenkendorf, R. and Mangold, M. (2014). Parameter Identification for Ordinary and Delay Differential Equations by Using Flat Inputs. Theoretical Foundations of Chemical Engineering, 48(5), 594–607. Schenkendorf, R., Reichl, U., and Mangold, M. (2012). Parameter Identification of Time-Delay Systems: A Flatness Based Approach. IFAC Proceedings Volumes, 45(2), 165–170. Verheijen, P.J. (2003). Model Selection: An Overview of Practices in Chemical Engineering. Computer Aided Chemical Engineering, 16, 85–104. Victor, S., Melchior, P., L´evine, J., and Oustaloup, A. (2014). Flat Output Computation for Fractional Linear Systems: Application to a Thermal System. IFAC Proceedings Volumes, 47(3), 2891–2896. Voit, E.O., Martens, H.A., and Omholt, S.W. (2015). 150 Years of the Mass Action Law. PLoS Computational Biology, 11(1). Wagner, M.O. and Meurer, T. (2018). Flatness-based feedforward control design and two-degrees-of-freedom tracking control for semilinear plug flow reactors. Journal of Process Control, 64, 132–140. Waldherr, S. and Zeitz, M. (2008). Conditions for the Existence of a Flat Input. International Journal of Control, 81(3), 439–443. Waldherr, S. and Zeitz, M. (2010). Flat Inputs in the MIMO Case. IFAC Proceedings Volumes, 43(14), 695– 700. Wey, T. (2002). Nichtlineare Regelungssysteme. Vieweg+Teubner Verlag, Wiesbaden.