Calculating D-optimal designs for compartmental models with a Michaelis–Menten elimination rate

Journal of Process Control 83 (2019) 88–101

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

Calculating D-optimal designs for compartmental models with a Michaelis–Menten elimination rate Belmiro P.M. Duarte a,c,∗ , Anthony C. Atkinson b , José F.O. Granjo c , Nuno M.C. Oliveira c a b c

ISEC, Department of Chemical & Biological Engineering, Polytechnic Institute of Coimbra, Rua Pedro Nunes, 3030-199 Coimbra, Portugal Department of Statistics, London School of Economics, London WC2A 2AE, United Kingdom CIEPQPF, Department of Chemical Engineering, University of Coimbra, Rua Sílvio Lima, Polo II, 3030-790 Coimbra, Portugal

a r t i c l e

i n f o

Article history: Received 8 June 2019 Received in revised form 3 September 2019 Accepted 3 September 2019 Keywords: Model-based optimal design of experiments Compartmental model Michaelis–Menten elimination rate Dynamic experiments Dynamic optimization

a b s t r a c t The optimal design of experiments for the parameters of the differential equations arising in many kinetic and pharmacodynamic models requires the numerical maximization of a function of the numerical solutions of sets of differential equations, the solutions depending on the experimental design. This problem in the design of experiments is solved as a dynamic optimization problem, using a simultaneous approach. The time horizon of the experiments is discretized in finite elements and orthogonal collocation on finite elements is used to parameterize the solution. The solution in each finite element is described by cubic Lagrange polynomials and the control action is represented by piecewise constant polynomials. We find D-optimum designs for two- and three-compartment models with Michaelis–Menten elimination rate kinetics. We consider three different design problems including both dynamic and static experiments. Single- and multi-experiment designs are addressed. Importantly, our algorithm allows the calculation of designs with constraints not only on the design region but also on the concentrations of the components in the compartments and on the rates of change of the design variables and of the resulting concentrations. © 2019 Elsevier Ltd. All rights reserved.

1. Motivation This paper addresses the D-optimal design of experiments for compartmental models which find application in many medical areas, in industrial processes and in systems theory. Here, we consider applications of interest in pharmacokinetics, when the compartments represent different sections of the body. The models consist of sets of differential equations describing the dynamics of the concentrations of a drug in the compartments over time. The unknown parameters of the models need to be estimated in order to determine the most suitable treatment regime, including dose levels. Because the economic resources available for experimental work are limited, optimally designed experiments are sought so that the information obtained is maximal. The methods of optimal experimental design are of particular importance in these models, which are generally nonlinear in inputs and parameters.

∗ Corresponding author at: Department of Chemical and Biological Engineering, ISEC, Polytechnic Institute of Coimbra, Rua Pedro Nunes, 3030-199 Coimbra, Portugal. E-mail addresses: [email protected] (B.P.M. Duarte), [email protected] (A.C. Atkinson), [email protected] (J.F.O. Granjo), [email protected] (N.M.C. Oliveira). https://doi.org/10.1016/j.jprocont.2019.09.001 0959-1524/© 2019 Elsevier Ltd. All rights reserved.

Compartmental models represent sets of interlinked wellmixed, homogeneous, cells within a system through which entities (e.g., material, energy) flow. They greatly simplify the mathematical description of complex systems in pharmacokinetics and epidemiology. In pharmacokinetics, they correspond to sets of ordinary differential equations (ODEs) describing the concentration dynamics of a drug within different sections of the body over time. Compartmental models are also useful in pharmacodynamic studies to establish the relationship at a site of action between the drug concentration and the resulting effect, the variation of the effect over time, the intensity of therapeutic and adverse effects. Compartment 1, in the pharmacodynamics conceptual framework, is where the drug is administered and the following ones, connected to the first, are where the drug effects occur [1]. The complexity of the ODEs modeling the drug concentration dynamics in compartments depends on the administration policy as well as on the elimination kinetic rate (clearance), which models the drug release from Compartment 1. While analytical solutions exist when the system of ODEs is linear, generally they are not available for nonlinear ODEs, requiring numerical methods to be used [2]. This is the case when clearance is modeled assuming a Michaelis–Menten elimination rate [3].

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

The optimal design of experiments (DOE) requires the numerical maximization of an objective function resulting from the solution of the ODEs. Atkinson and Bogacka [4] calculated D-optimal DOE to estimate the parameters of kinetic models, aiming to find the optimal set of sampling times from a chemical reaction apparatus. Atkinson and Fedorov [5] introduced T-optimum experimental designs for discriminating between models. The generalization by ´ Ucinski and Bogacka [6] to designs for discriminating between chemical kinetic models extends the solution from the establishment of time points to include the calculation of an optimal temperature profile. A simulation-based method was adopted by Franceschini and Macchietto [7] to determine not only the sampling instants but also the initial conditions and the optimal policy of control actions, such as the dose as a function of time and, in chemical experiments, the temperature profile. Similarly, Galvanin et al. [8] used online model-based redesign of experiments and backoffbased optimal designs to prescribe optimal dynamic experiments for parameterizing a type 1 diabetes model. Other methods for model-based optimal design of experiments (MbODE) applied to dynamic models treat the problem as an optimal control problem and use dynamic optimization methods to ´ provide numerical answers, e.g., Ucinski [9, Sec. 4.2] and Pronzato [10]. Approaches handling the dynamic optimization problems include the sequential [11,12], multiple-shooting [13], and the simultaneous [14] methods. Sequential methods use a control vector parametrization technique, where the input trajectory is parameterized over the time interval [15,16]. Multiple-shooting techniques divide the time interval into smaller intervals, solving an initial value problem in each of the sub-intervals while imposing constraints that ensure the properties of the solution in the original time horizon [17–19]. In simultaneous approaches, the time horizon is discretized into finite elements on top of which are constructed algebraic approximations of the ODEs solution. Orthogonal collocation on finite elements (OCFE) is used for local approximation of the variables and Lagrange polynomials are employed to approximate the solution of each variable on each finite element. More details on OCFE can be found in Carey and Finlayson [20]. These approaches avoid the need for embedding a differential algebraic equation (DAE) solver in the optimizer and produce problems with sparsity and structure that can be exploited, improving the computational efficiency. These properties make simultaneous approach appropriate for the optimal design of experiments. But, from our knowledge, the only applications are Hoang et al. [21] who apply the approach to a semi-batch reactor and Duarte and Atkinson [22, Sec. 5] where they address temperature-dependent kinetic rates. This paper addresses the D-optimal design of experiments for compartmental models with a Michaelis–Menten elimination rate, handling both static and dynamic kinds of experiments. Static experiments are those where inputs are initially chosen and kept constant over time, while in dynamic experiments the inputs can vary during the experiment at a previously defined grid of time instants. In all scenarios the process is also sampled at a previously set grid of times, that may (or may not) coincide with the grid used to manipulate the inputs (modifying the dose) in dynamic experiments. The problems to be tackled are more complicated than those solved by, for example, Atkinson and Bogacka [4] who use the direct method of Valko and Vajda [23] to calculate parameter sensitivities from the numerical solution of sets of differential equations. These are in turn used to construct the Fisher information matrix (FIM) for finding optimal experimental designs. Here, in addition to the constraints that define the design region, constraints on control variables and on the concentrations of components and their gradients in the compartments can be included in the problem. Consequently, the dynamic optimization problem of maximizing the information from experiment(s) can have multiple local optima,

89

and a global optimizer must be used to guarantee the global optimal design. Another novelty in this work is the use of a “simultaneous” method to find D-optimal designs for compartmental models with Michaelis–Menten elimination rate. Matrix-based operations are employed to construct an explicit formulation of the optimality criterion, which is embedded in the design problem. We focus on D-optimality and use the Cholesky decomposition of the FIM for formulating the design problem, addressing different contexts where various kinds of action are possible. Additional points of innovation include (i) the generalization of the formulation from static to dynamic experiments; (ii) the generalization from single- to multi-experiment designs; and (iii) the generalization from unconstrained to constrained designs where both the control and the state variables can be bounded. The structure of the compartmental models studied herein is similar to that of the models of Galvanin and Bezzo [24] where they applied optimal DOE to in vitro bacterial activity and von Willebrand disease dynamics, and Telen et al. [25] where they use a multiple-shooting approach for finding optimal designs of experiments for nonlinear dynamic (bio)chemical systems. The framework here proposed is broad enough to be successfully applied to more complex problems including pharmacodynamic models of interest. The paper is organized as follows. Section 2 introduces the background and the notation used to formulate the problem as well as the fundamentals of nonlinear programming. Section 3 presents the mathematical programming formulation for finding the optimal designs for all the design setups addressed: 1. experiments where the only decision variable (experimental factor) is the input level of the drug, which is kept constant for each run; 2. experiments where the decision variables are the input level and the initial concentration in Compartment 1; and 3. experiments where the inputs vary during the experimental run. First, formulations are developed for single-experiment optimal designs, and later extended to combined multi-experiments setups. Section 4 applies these formulations to find D-optimal designs for two- and three-compartment models in the absence of constraints on temporal changes in the control variables. Section 5 demonstrates the effect of three kinds of constraints: 1. constraints on the magnitude of changes of the control variables; 2. constraints on changes in responses; and 3. constraints on responses themselves. Finally, Section 6 offers a summary of the results obtained and Appendix A provides technical details of the time discretization procedure. 2. Notation and background The nomenclature used in the representation of the models is introduced herein, and in Section 2.1 the experimental design problems sketched above are defined. Bold face lowercase letters represent vectors, bold face capital letters continuous domains, blackboard bold capital letters discrete domains and capital letters matrices. Finite sets containing  elements are compactly represented by [] ≡ {1, . . ., }. The transpose operation of a matrix is represented by “ ”. The concentration of active agent in compartment i, and its volume, are denoted by Ci and Vi , respectively. The mass transfer rate constant from compartment i to compartment j is ki,j , u1 (t) is the

90

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

stimulus applied in Compartment 1, where u1 (t) can vary with time and represents the policy of administration of the active principle. Practically, u1 (t) aggregates the information relative to dose and administration policy (i.e., the time dependency). m denotes the maximum rate achieved at saturating substrate concentration in Michaelis-Menten elimination rate models (from Compartment 1) and  m , the substrate concentration at which the elimination rate is m /2. Note that only C1 (t) is measured in all experiments. Consider a generic model with c compartments where the concentrations in the compartments, Ci (t), i ∈ [c], are process states, and the input u1 (t) is the control variable, varied over time in dynamic experiments, but kept fixed in static experiments. The values of u1 (t) are chosen to maximize the information obtained from the experiments. Let the vector of states be s(t) ≡ {Ci (t), i ∈ [c]}, the vector of process measurements, denoted as y(t), containing a single element, i.e., y(t) ≡ {y1 (t)} where y1 (t) = C1 (t) + (t) and (t) is the observational error, the vector of control variables be u(t) ≡ {u1 (t)}, and the vector of parameters  including all the parameters in the model that are to be estimated from the experiment(s), with s(t) ∈ S ⊂ Rns , y(t) ∈ Y ⊂ Rny and u(t) ∈ U ⊂ Rnu .  [ LO ,  UP ] is a comThe parameters  ∈ P ⊂ Rn where P = ⊗nj=1 j j pact set in the domain of the parameters,  j represents a local value of parameter j and jLO and jUP are the lower and upper values admissible for  j , respectively. Further, p is a specific vector of parameters and is used to represent the values of  used in designing the experiment. Here, ns stands for the number of states (c), ny for the number of concentrations measured (i.e., one), nu for the number of control variables (one), n for the number of parameters, and S, Y, U and P are the compact domains of states, measurements, controls, and parameters. The two-compartment model (2-COMP) is first presented and then generalized to the three-compartment one (3-COMP). A representation of two- and three-compartment models is shown in Fig. 1. The 2-COMP model is as follows (1a)

k1,2 V1 dC2 = C1 − k2,1 C2 , V2 dt

(1b)

C1 (0) = C1in ,

(1c)

y1 (ti ) = C1 (ti ) + (ti ),

(1d) g

ti ∈ t ,

(1e)

where u1 (t) is the variation in drug concentration fed to Compartment 1 per unit of time with C1in and C2in the drug concentration in Compartments 1 and 2 when the experiments start. Eqs. (1a) and (1b) represent the dynamics of drug concentration in compartments 1 and 2; Eq. (1c) is the initial condition for C1 and Eq. (1d) for C2 . Here, tg is the set of time instants ti , i ∈ [nt ] at which the time domain is discretized and measurements are available. The number of time instants is nt and ti ∈ [0, H], H being the horizon of the experiment(s). We consider the measurements are discretely taken at instants tg and this grid is simultaneously used for discretizing the time horizon of the experiments by choosing the size of the finite elements. Thus, Eq. (1e) is the measurement equation where the errors (t) are assumed independent and identically distributed with mean zero and constant variance. Next, the 3-COMP model is given by:

k2,1 V2 k3,1 V3 C2 − k1,3 C1 + C3 + +u1 (t), V1 V1

k1,3 V1 dC3 = C1 − k3,1 C3 , V3 dt

(2c)

C1 (0) = C1in ,

(2d)

C2 (0) = C2in ,

(2e)

C3 (0) = C2in ,

(2f) g

y1 (ti ) = C1 (ti ) + (ti ),

ti ∈ t .

(2g)

C3in

is the drug concentration in Compartments 3 when the Here, experiments start. For compactness, 2-COMP model (Eq. (1)) is designated as M2 (y, s, u, ) and the 3-COMP model (Eq. (2)) as M3 (y, s, u, ). The discretized forms in time domain of 2-COMP and 3-COMP models are expressed as Md2 (y, s, u, ) and Md3 (y, s, u, ), respectively. Appendix A details the discretization procedure. The systems of ODEs describing the behavior of the concentration in compartments need to be integrated numerically together with the equations for the parameter sensitivities. The continuous method for computing the sensitivity equations is used, since is known to be accurate and not too computationally demanding for small problems [26]. Calculation of parameter sensitivities via the continuous method can be found in the analysis of kinetic systems; see Saltelli et al. [27] or Maly and Petzold [28]. The discretization of the continuous sensitivity equations uses the same scheme as we apply to the ODEs of the compartmental models. Let the sensitivity of the state variable si with respect to parameter  j be denoted by  i,j . The sensitivities of the states with respect to the parameters are obtained by solving the additional ODEs di,j dt

=

ns  ∂f (s, u, ) i

∂sk

i,j (0) = 0,

k,j +

∂fi (s, u, ) , ∂j

i ∈ [ns ], j ∈ [n ] (3a)

i ∈ [ns ], j ∈ [n ],

(3b)

generically denoted by S2 (, y, s, u, ), together with the model equations. Here, fi (s, u, ) is the right hand side function of the ODE for state si , i.e., for 2-COMP models f1 (s, u, ) = − V (m+C ) C1 − k

V

k

V

1

m

1

2 1 C2 + u1 (t) and f2 (s, u, ) = 1,2 C1 − k2,1 C2 . Equivak1,2 C1 + 2,1 V1 V2 lent relations can be established for the 3-COMP models. The sensitivity of the general measurement variable yk with respect to parameter  j is represented by ϑk,j , i.e.,

ϑk,j =

ns  ∂gk (s, u, )

∂si

i=1

i,j +

∂gk (s, u, ) , ∂j

k ∈ [ny ], j ∈ [n ], (4)

where gk (s, u, ) is the measurement equation for response yk . Since for both models there is only a single measurement, g1 (s, u, ) = C1 (t). The set of parameters to be estimated from the experiment(s) consists of both parameters of the Michaelis–Menten rate of elimination and the rate constants ki,j . When 2-COMP models are considered,  ≡ {m ,  m , k1,2 , k2,1 }. Here, the volumes Vi , i ∈ [c] are assumed known a priori, and therefore not included in . Without loss of generalization the volumes can also be considered as parameters to find and included in . Applying (3), the sensitivity equations for 2-COMP model yield





d1,1 k2,1 V2 m m =− + k1,2 1,1 + 2,1 V1 dt V1 (m + C1 )2

m dC1 =− C1 − k1,2 C1 dt V1 (m + C1 ) +

(2b)

k=1

k2,1 V2 dC1 m =− C1 − k1,2 C1 + C2 + u1 (t), V1 dt V1 (m + C1 )

C2 (0) = C2in ,

k1,2 V1 dC2 C1 − k2,1 C2 , = V2 dt

(2a)

−

C1 , V1 (m + C1 )

(5a)

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

91

Fig. 1. Systematic representation of (a) 2-COMP model and (b) 3-COMP model.





k2,1 V2 d1,2 m m + k1,2 1,2 + 2,2 =− V1 dt V1 (m + C1 )2 +

m C1 V1 (m + C1 )2

,



(5b)



d1,3 k2,1 V2 m m =− + k1,2 1,3 + 2,3 − C1 , 2 V1 dt V1 (m + C1 )

(5c)

d1,4 k2,1 V2 V2 m m =− + k1,2 1,4 + 2,4 + C2 , V1 V1 dt V1 (m + C1 )2

(5d)

d2,1 k1,2 V1 1,1 − k2,1 2,1 , = V2 dt

(5e)

d2,2 k1,2 V1 1,2 − k2,1 2,2 , = V2 dt

(5f)

d2,3 k1,2 V1 V1 1,3 − k2,1 2,3 + C1 , = V2 V2 dt

(5g)

d2,4 k1,2 V1 1,4 − k2,1 2,4 − C2 , = V2 dt

(5h)





i,j (0) = 0,

i ∈ [ns ], j ∈ [n ],

(5i)

ϑ1,j = 1,j ,

j ∈ [n ],

(5j)

where (5i) represents the initial conditions for the sensitivities and (5j) the sensitivities of the measurement which are, in turn, used for constructing the FIM. The sensitivity equations for 3-COMP models are similarly constructed and their development is omitted. The ODEs system (5) is compactly represented by S2 (, y, s, u, ) and the sensitivities for 3-COMP model are denoted as S3 (, y, s, u, ). Their discretized representations are Sd2 (, y, s, u, ) and Sd3 (, y, s, u, ), respectively. The extension of the approach to generic models formalized by DAEs is straightforward; in every case the sensitivity computation requires solving an additional set of ODEs to determine . To avoid any adverse numerical effects of the magnitudes of the sensitivities we apply a scaling procedure to ϑk,j , k ∈ [ny ], j ∈ [n ] and in the following mathematical programming formulations work with the scaled sensitivities  [21], i.e.,

k,j = pj ϑk,j ,

k ∈ [ny ], j ∈ [n ].

(6)

The solution of 2-COMP models is carried out simultaneously with the solution of the sensitivity equations Sd2 (, y, s, u, ) and is detailed in Appendix A. Noticeably, the numerical solution of the ODEs system using OCFE is enclosed in an optimization problem of nonlinear programming (NLP) type that maximizes a measure of information for the experiment defined by a specific convex function of the FIM. Let experimental designs be denoted by  and formed by a vector of actions. When errors are identically and normally distributed

with independent components, the volume of the confidence region for p is inversely proportional to det [M1/2 (|u1 , tg , p)] where the FIM is denoted by M(|u1 , tg , p) and depends on the input variables and the grid of discretized time instants, tg . Consequently, minimizing the determinant of the inverse of the FIM by choice of u1 in each discretized time interval leads to the most accurate estimates for the parameters. Other alphabetic criteria optimize the information content in different ways [29,30] but in practice the D-optimality criterion is the most used. In subsequent sections we consider equally spaced time grids but other schemes can be used. In this study, additive errors of constant variance in the measurements are assumed. Other distributional assumptions about the observational errors lead to other forms for the FIM [31]. 2.1. Optimal design problems Let us consider single-experiment setups where the parameters are to be fit from a single experiment, and the decision variables of the design problem are chosen so that the information obtained is maximum. The three kinds of D-optimum design problems to be addressed are as follows: 1. SE – Static Experiments – here the experimental or decision variable is the initial value for the drug administration rate, u1 (0), which is kept fixed during one experimental run, i.e., u1 (0) = u1 (t), t ∈ [0, H]. Consequently, the design problem consists of finding the optimal value of u1 (0) so that det [M(|u1 (0), tg , p)] is maximized. 2. SE&I – Static Experiments with, in addition, an initially chosen concentration in Compartment 1 – the decision variables are the initial drug administration rate, u1 (0), and the initial concentration in Compartment 1, i.e., C1 (0). Now we have to find the optimal values of u1 (0) = u1 (t), t ∈ [0, H] and C1 (0), so that det [M(|u1 (0), C1 (0), tg , p)] is maximized. 3. DE – Dynamic Experiments with chosen initial concentration in Compartment 1 – the decision variables are the initial drug administration rate, u1 (0) and the optimal profile of u1 (ti ), ti > 0, ti ∈ tg , for a discrete grid of times at which the dose may be optimally changed. We also choose the value of C1 (0). Now, det [M(|u1 (ti ), C1 (0), tg , p)] is maximized. Let vector z aggregate the decision variables, then for the SE design problem z ≡ {u1 (0)}, for SE&I design problem z ≡ {u1 (0), C1in }, and for DE, z ≡ {u1 (t), t ∈ t g , C1in }. All design problems have the drug concentrations in compartments and parameter sensitivities as decision variables but this dependency is omitted for simplicity. For the three kinds of design we have =





u1 (0) ,



=

u1 (0) C1in





,

and

=

u1 (ti )

C1in

ti

t0



,

ti ∈ t g

92

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

for problems 1, 2 and 3 above. Optimal designs are represented as * and t0 = 0 in all the experiments addressed. To compare the information content obtained from two differref , where the latter one is the reference, ent designs, say  D and D we use the D-optimality efficiency [29,Chap. 11]



EffD =

1/n

det[M(D , p)] ref , p)] det[M(D

.

(7)

NLP formulation used for determining D-optimal designs in the following sections reads

log[det{M(|z, t g , p)}]

(8a) (8b)

Sd2 (, y, s, u, p).

(8c)

Let each element of M(|z, tg , p) be designated as mi,j , i, j ∈ [n ]. Supposing y1 (t) as the only response measured, each vector of parameter sensitivities at time tl ∈ tg is (tl ) ∈ R1×n , l ∈ [nt ] and each element of (tl ) is i (tl ), l ∈ [nt ], i ∈ [n ]. Then, each element of the FIM is given by 1

i (tl ) j (tl ). nt nt

mi,j =

s.t. Model :

(9)

mi,j =

(12b)

Sd2 (, y, s, u, p) and (6),

1

1,i (tl ) 1,j (tl ), nt

(12c)

mi,j =



ai,l aj,l ,

ai,i ≥ p , ai,j = 0, mi,i ≥

a2i,j

(11)

n

 Equivalently, log[det{M(, z, t g , p)}] = 2 i=1 log(ai,i ), and g maximizing det {M(z|t , p)} is equivalent to maximizing the sum of the logarithms of the diagonal elements of A(z|tg , p). The

(12g)

i, j ∈ [n ],

(12h)

where p is a small positive constant to ensure the semipositiveness of the FIM. In all examples presented in the following sections we set p = 10−8 . Eq. (12d) is to construct the FIM, (12e) represents the Cholesky decomposition, (12f) is to guarantee the positiveness of the diagonal elements of A(z|tg , p), and (12g) ensures that A(z|tg , p) is upper diagonal. Eq. (12h) is a numerical stability condition imposed on the Cholesky factorization of positive semidefinite matrices [34, Theorem 4.2.8].. For the 3-COMP model, Md2 (y, s, u, p) in (12b) is replaced by Md3 (y, s, u, p), and the set of equations Sd2 (, y, s, u, p) in (12c) by Sd3 (, y, s, u, p). The strategy used to compute the determinant of the FIM, which is based on Cholesky factorization, avoids the explicit use of the det(•) function which is not available in the GAMS software we used for solving the optimization problems. This strategy is not required for software including the determinant function when (12a) is replaced by det {M(z|tg , p)}. Further, Eqs. (12e)–(12h) are not required, being implemented at a lower level of computation and assured by the function itself. Now, we generalize the formulation (12) for multi-experiments. We recall that ne designates the number of experiments to be carried out which is previously set. The constraint (12d) modeling the FIM construction is replaced by 1 

i,k (tl ) j,k (tl ), nt ne

 = (u1,k (0),

a2i,i .

(12e) (12f)

i, j ∈ [n ], i < j,

M(|z, t g , p) = A (z|t g , p)A(z|t g , p),

n 

i, j ∈ [n ], i ≤ j,

i ∈ [n ],

ne

where A(z|tg , p) is an upper triangular matrix with positive diagonal elements [33]. Afterwards, A(z|tg , p) is used to calculate the determinant of M(|z, tg , p). Let the elements of A(z|tg , p) be ai,j , i, j ∈ [n ], where ai,j = 0, ∀ j ≤ i + 1, i, j ∈ [n ]. The determinant of the FIM can be calculated from the diagonal elements of A(z|tg , p) using the equality

(12d)

l=1

mi,j =

(10)

i, j ∈ [n ],

l=1

This can be generalized for a different measurement, say of C2 , as well as to multiple measurement systems, where the covariance matrix of the responses is required [32]. Maximizing log [det {M(|z, tg , p)}] is equivalent to maximizing [det {M(|z, tg , p)}]. The semi-definite positiveness of the FIM can be exploited by applying the Cholesky decomposition

i=1

Md2 (y, s, u, p),

n

l=1

det{M(|z, t g , p)} =

(12a)

nt

s.t. Model Md2 (y, s, u, p), Equations

log(ai,i )

i=1

Sensitivity equations :

In this section we describe the numerical procedure for finding optimal experimental designs for estimation of the parameters  = (m ,  m , k1,2 , k2,1 ) in 2-COMP models. The extension to 3-COMP models, where the set of parameters is  = (m ,  m , k1,2 , k2,1 , k1,3 , k3,1 ) is straightforward, as it is to any other compartment model. First, we consider single-experiment designs and provide formulations for this setup. Next, the extension to multi-experiments designs for SE and SE&I setups is addressed as it results from generalizing the previously introduced problems. We notice that the goal for multi-experiment design problems is choosing the optimal profiles of the control variable for a combined set of ne experiments with the aim of maximizing the overall information obtained. The main motivation for analyzing multi-experiment setups is that they assure higher efficiencies than a single experiment, obviously with additional costs, and are advantageous over single-experiment designs when the control factors that can be changed are the level of the control variable at the beginning and the initial concentration, as for the SE and SE&I design problems. The D-optimal design problem on U × S is max

2

z ∈ U,s ∈ S,y ∈ Y

3. Finding D-optimal designs

z ∈ U,s ∈ S,y ∈ Y

n 

max

nt

i, j ∈ [n ].

(13)

k=1 l=1

Here, i,k (tl ) is the sensitivity with respect to the ith parameter at time tl ∈ tg for kth experiment, and the optimal design is represented by k ∈ [ne ]),

where u1,k (0) is the input value used in kth experiment. The model Md2 (y, s, u, p) and the sensitivity equations Sd2 (, y, s, u, p) are solved simultaneously for ne experiments using the same discretization grid, tg . In each experiment the control variable profile and/or the initial condition are chosen so that the combined amount of information is maximized. In the absence of any constraints on the problem, other than those due to the design region, the maximization of the determinant of the FIM as a function of the approximate design measure  is a problem in convex optimization. However, it becomes a rather challenging NLP problem due to the need to compute the determinant via Cholesky factorization. This produces multinomial terms of different orders of magnitude which often lead to ill-conditioned

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101 Table 1 Parameters for constructing locally D-optimal designs. Parameter

Value

Parameter

Value

Parameter

Value

k1,2 k3,1 V1

1.0 h−1 0.8 h−1 1.0 L

k2,1 m V2

1.0 h−1 0.2 h−1 3.0 L

k1,3 m V3

0.3 h−1 0.01 g/L 1.5 L

Table 2 Variable bounds and other information required by the formulation. Time horizon Bounds for input variable Bounds for state variable Finite elements Local discretization points tg

H = 10 h u1 (t) ∈ [0.05, 0.5] g/(L h) Ci (t) ∈ [0, 1.5] g/L, i ∈ [c] E = 20 (of equal size) K = 4 (one coinciding with the bound) formed by all the finite elements bounds and discretization points

93

All the following numerical examples consider the absence of drug in the second and third compartments at the beginning of the experiments, i.e., Ciin = 0, i ∈ {2, . . ., c}. This is typically the most realistic scenario, but other values could be used. Herein, we consider an acquisition frequency of 2 measurements/h, a rate currently ensured by common instrumentation [39]. The results obtained for 2-COMP and 3-COMP models and for the three design scenarios (SE, SE&I and DE) described in Section 2.1 are now presented. We recall that the set of parameters to be estimated for the 2-COMP model is  ≡ {m ,  m , k1,2 , k2,1 } while for the 3-COMP model two additional parameters are included, i.e.,  ≡ {m ,  m , k1,2 , k2,1 , k1,3 , k3,1 }. We notice that the optima reported for all the designs is the objective in problem (12), i.e., the value of log [det {M(|z, tg , p)}] used for computing the D-optimal efficiency and the determinants ratio. 4.1. SE: optimal design for Static Experiments

FIMs. Especially important is the presence of constraints on the values of the concentrations and of the rate of change of the concentrations and the control variables, which we consider in Section 5. Since the optimal design problem (12) may have multiple local optima we require a global solver. To determine the optimal design we coded the problem in GAMS [35] and used a multistart heuristic algorithm-based solver, OQNLP, to find the global optimum. The algorithm calls an NLP solver from multiple starting points, keeps all the feasible solutions found, and picks the best as the optimal solution of the problem [36]. The starting points are computed with a random sampling driver that uses independent normal random variables for each decision variable. Practically, OQNLP does not guarantee that the final solution is a global optimum, but it has been successfully tested on a large set of global optimization problems. To build the initial sampling points the variables need to be bounded, which is what we have since the design space and the region of plausible values are all compact by assumption. The NLP solver called by OQNLP is CONOPT, which in turn uses the generalized reduced gradient (GRG) algorithm [37]. The maximum number of starting points allowed is set to 3000 and the procedure terminates when 100 consecutive NLP solver calls result in an improvement less than a tolerance of 10−4 . The absolute and relative tolerances of the solvers were all set equal to 10−8 and 10−7 , respectively. All computations in Section 4 used an Intel Core i7 machine running a 64 bits Windows 10 operating system with a 2.80 GHz processor.

4. Unconstrained D-optimal designs This section presents the D-optimal designs calculated when the only constraints on the decision variables (inputs and states) are physical bounds on U × S. In Section 5 we find designs for dynamic experiments when there are, in addition, constraints on the values of u1 (t) and C1 (t). The information required for the scenarios listed in Section 2.1 is in two parts. Table 1 contains the parameters used for constructing locally D-optimal designs. The values of the parameters used in the numerical experiments were chosen to give ratios of ki,j rates similar to those found in several studies available in the literature. Table 2 presents other data required by the algorithm, namely the bounds used for the decision variables and the discretization grid. Here, we adopt an equidistributed grid of points (equal size finite elements) since previous simulation studies indicated that the solution activity estimated from the second order derivatives of the state variables is of the same order of magnitude over the horizon of the experiments [38], i.e., the local error of the ODEs solution is similar all over the experiment(s) duration.

The decision variable is u1 (0), the initial condition for concentration in Compartment 1 is C1 (0) = 0.0 g/L, and the complete experiment generates 21 measurements coinciding with the bounds of the time intervals (i.e., finite elements). Fig. 2(a) and (b) displays the dynamics of the concentrations over time for the 2COMP and 3-COMP models. We notice that the profiles are similar; all concentrations have an increasing behavior as a result of the high prescribed value of u(0). The increasing behavior for the concentrations was expected given the initial condition for C1 (0). The optimal value of u(0) is limited by the upper bound on C1 at the end of the experiment for the 2-COMP model (see Fig. 2(a)) and by the upper bound of u(0) itself for the 3-COMP model (see the second line of Table 3). Optimal plans for two- and three-experiment designs were obtained to analyze the impact of the number of experiments on the efficiency. The D-optimal efficiency of the designs shown in Table 3 are computed having the three-experiment design for reference ref ) and using Eq. (7). The optimal design for the four-experiment (D plan was also found but is not included since its efficiency is similar to that of the optimal three-experiment plan. The efficiency of the single experiment plan is so low that these designs have no value. Practically, to extract valuable information in this scenario we need more than one experiment. The results show that the efficiency increases as the number of experiments of the plan increases when ne ≤ 3 for both the 2-COMP and 3-COMP models which include 4 and 6 parameters, respectively. The upper bound of u1 (0) is a support point in all calculated designs, while the other support points change. Fig. 3(a) and (b) displays the dynamics of the concentrations in the compartments for 2-COMP and 3-COMP models for a two-experiment design. Numerical tests revealed that, in these examples, the effect of scaling the sensitivity coefficients via Eq. (6) is negligible, because the variables ϑk,j are of a similar order of magnitude. Nonetheless, for general problems the scaling might minimize numerical problems in the FIM construction and subsequent Cholesky decomposition. 4.2. SE&I: optimal design for Static Experiments with choice of initial concentration The decision variables are u1 (0) and the initial concentration in Compartment 1, C1 (0). The optimal designs for multi-experiment setups are represented by



=

u1,k (0) C1,k (0)



,

k ∈ [ne ],

94

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

Fig. 2. SE: concentration dynamics resulting from optimal single-experiments for: (a) two-compartment model; (b) three-compartment model.

Table 3 SE: Optimal design for Static Experiments plan. Model

ne

Optimal design

2-COMP 3-COMP

1 1

2-COMP

2

(0.4810)  = (0.5000) 

3-COMP

2

2-COMP

3

=

3-COMP

3

=

a

=

=

 0.2303

 0.4810 

0.2525 0.5000

  0.1793

0.2809 0.4810

0.1997 0.3090 0.5000

 

Optimum

EffD a

CPU time (s)

−10.7552 −18.8060

0.2542 0.1453

349.92 432.12

−5.5659

0.9301

435.78

−10.1517

0.6149

613.32

−5.2759

1.0000

673.42

−7.2338

1.0000

859.85

Computed from Eq. (7) with * being the 3-experiment designs.

Fig. 3. SE: concentration dynamics resulting from optimal two-experiments obtained for: (a) two-compartment model; (b) three-compartment model.

where the first line contains the inputs in each experiment, and the second the initial condition. First, we address the optimal designs for one-experiment plans. Fig. 4(a) and (b) presents the evolution of the concentrations for optimal plans with a single experiment for the 2-COMP and 3-COMP models, respectively. The optimal profiles show differences relative to those obtained for the SE designs. Here, the initial condition is close to its upper value and C1 (t) decreases with time. Contrarily, the control variable is at its lower bound. Table 4 presents the optimal designs generated and the efficiencies obtained with the three-experiment design as reference.

Fig. 5 shows the behavior of the concentrations in the compartments for two-experiment plans, and we notice that two trend curves appear for C1 (t). For C1 (0) = 1.5 g/L, C1 (t) initially decreases and then increases again to a final value of 1.5 g/L at the end of the run; the value of u1 (t), well below the boundary of the experimental region ensures that Ci (t) ≤ 1.5 g/L for all t. On the other hand, for C1 (0) close to zero, C1 (t) increases steadily, as do the two curves for C2 (t) and C3 (t), see Fig. 5(a) and (b). Like the SE design, here the efficiency of the designs also increase as the number of experiments increases for ne ≤ 3. However, unlike the SE design, the extra information from the choice of C1 (0) makes the single experiment design

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

95

Fig. 4. SE&I: concentration dynamics resulting from optimal single-experiments obtained for: (a) two-compartment model; (b) three-compartment model. Table 4 SE&I: optimal design for Static Experiment plans with choice of initial concentration. Model

ne

2-COMP

1

3-COMP

1

2-COMP

2

3-COMP

2

2-COMP

3

=

3-COMP

3

=

a

Optimal design

=

1.4826 

=

= =

0.0500



0.0500

1.4598  0.2134 0.3412

0.0039

1.5000  0.1138 0.3924

0.6688 1.5000  0.0845 0.1942 0.3412

1.5000

0.0093 1.5000  0.0826 0.2155 0.3924 1.5000

0.0043

1.5000

Optimum

EffD a

CPU time(s)

−4.5908

0.6739

310.23

−11.1368

0.2961

328.28

−3.5784

0.8059

397.34

−6.9713

0.5929

475.44

−3.0124

1.0000

534.28

−3.8348

1.0000

978.22

Computed from Eq. (7) with * being the 3-experiment designs.

Fig. 5. SE&I: concentration dynamics resulting from optimal two-experiments obtained for: (a) two-compartment model; (b) three-compartment model.

reasonably efficient for the 2-COMP model (67.39%), although not for 3-COMP. 4.3. DE: optimal design for dynamic experiments Here, we consider the optimal design for a dynamic experiment where the control action, the input dose u1 (t), can be manipulated,

and C1in (equal to C1 (0)) chosen so that the information obtained from a single experiment is maximized. We recall that the action is described by piecewise constant polynomials in each finite element and the state variables, as well as the sensitivities, by cubic Lagrange polynomials. To represent the control variable u1 we assume that it is constant in each finite element and its value results from the optimization procedure. The changes only occur at the bounds of

96

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

Fig. 6. DE: results for 2-COMP model: (a) optimal profile of dose u1 (t); (b) dynamics of the concentrations in both compartments.

the finite elements (time intervals previously set) and are described by steps. Specifically, we have

u1 (t) =

⎧ u1 (0) ⎪ ⎪ ⎪ ⎪ ⎪ u1 (t1 ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ··· ⎪ ⎪ ⎪ ⎨ u1 (t )

if t ∈ [0, h1 [, if t ∈ [h1 , h1 + h2 [, ··· if t ∈



i 

h,

i+1 

 h

,

i i i ⎪ ⎪ i=1 i=1 ⎪ ⎪ ⎪ ⎪ ··· ··· ⎪ ⎪ NE−1 NE  ⎪ ⎪   ⎪ ⎪ ⎪ hi , hi . ⎪ ⎩ u1 (tNE−1 ) if t ∈ i=1

(14)

i=1

To explicitly formalize the discontinuity at the nodes and the Heaviside nature of the control variable we use open intervals at the end of the finite elements and closed intervals at the beginning. The exception is tUP where Eq. (A.5) holds. The assumption that the control variables can be manipulated instantaneously may be challenging in practical applications, but currently the actuation time is in most of the cases negligible compared to dominant process dynamics time constants, and the optimal policy prescribed with our formulation is acceptably close to the optimum. If finer detail is required for input representation, higher degree polynomials can be used, which have the drawback of increasing the size of the optimization problem, since more variables must be optimized per finite element. Nonetheless, the extension to higher polynomial representations of the formulation proposed here, which may allow avoiding instantaneous actions, is straightforward. Our strategy assumes the optimal profile of actions is implemented, and if distributional limitations and/or adverse effects on other variables occur, they can be included in the model and would lead to additional constraints on the optimal design problem. Section 5 demonstrates the application of the proposed approach to setups where constraints are explicitly considered. In this section we limit our analysis to single-experiment scenarios since the quantity of resulting measurements and the excitation induced by the control variable are generally rich enough for fully parametrization of the model. Multi-experiment setups can be addressed similarly, as in Sections 4.1 and 4.2. Fig. 6(a) shows the action profile for the 2-COMP model where u1 (t) twice visits the upper bound of 0.5 g/(L h) as the experiment evolves. The dynamics of the concentrations are shown in Fig. 6(b) and reveal an initial increase of C1 followed by decrease and a final

increase. The results for the 3-COMP model are in Fig. 7(a) and (b). Here, u1 (t) only visits the upper bound once. Numerical results for the optimal designs are in Table 5 where the first line of  is for u1 (t) and the second for intervals where it is constant. The last column is for optimal C1in (in the first line) and t0 (in the second line), which is 0.0 h in all experiments. It is well-known that D-optimal designs for linear models choose experimental points at the extremes of the design region. The same is found [40] for D-optimal designs for linear control systems, resulting in persistent excitation. Here, the optimal profiles of actions likewise require extreme variations of the input, going suddenly from lower to upper bounds or vice-versa. Note that, although our models are not linear, some parameters do occur linearly in the kinetic rate equations of Section 2. In cases when the magnitude of the changes in control or state variables is limited, constraints in the optimal design problem can be included, as is shown in Section 5. The advantage of the NLP formulation proposed is that it allows a straightforward inclusion of the constraints in the design problem. 5. DEC: optimal designs for dynamic experiments with constraints Formulation (12) is used to find constrained D-optimal dynamic designs for the 3-COMP model (2). First, we address the design problem where the changes in u1 (t) between two consecutive discrete times of the grid tg are constrained to a pre-specified interval [− u1 , + u1 ], with u1 the maximum change allowed for control variable u1 . The design problem is solved including the constraints u(ti ) − u(ti−1 ) ≤ + u1 ,

ti ∈ t g , i ≥ 1,

(15a)

u(ti ) − u(ti−1 ) ≥ − u1 ,

ti ∈ t g , i ≥ 1

(15b)

in program (12). Let u1 = 0.20 g/(L h) and call the optimal design obtained  constr,u . The optimal profile of actions of  constr,u is in Fig. 8(a). The difference with the unconstrained design in Fig. 7(a) is small. Specifically, three sequential steps of the control variable are required to cover the complete range of C1 , rather than two. The behavior of the concentrations in all compartments is only slightly affected as can be observed on careful comparison of Fig. 8(b) and Fig. 7(b). The first line of Table 6 shows the constrained optimal design  constr,u . Its efficiency computed with equation (7) and using the equivalent unconstrained design for reference (in the second line

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

97

Fig. 7. DE: results for 3-COMP model: (a) optimal profile of dose u1 (t); (b) dynamics of the concentrations in the three compartments. Table 5 DE: optimal design for dynamic experiment. Model Two-compartment Three-compartment

Optimal design

=

0.1924,

0.5000,

0.0500,

0.5000

0.0100



[0.0, 1.0[, [1.0, 4.5[, [4.5, 8.5[, [8.5, 10.0[

0.0 0.0500, 0.1577, 0.5000 1.008 = [0.0, 6.5[, [6.5, 7.0[, [7.0, 10.0[ 0.0

Optimum

CPU time (s)

−5.5725

238.39

−9.2566

353.54

The first line of  is for u1 (t) and the second for time intervals where it is constant. The last column is for C1in and t0 .

Fig. 8. DEC: dynamic optimal design for 3-COMP model with constraints (15) on input variation: (a) optimal profile of actions; (b) dynamics of the concentrations in the compartments.

of Table 5) is 99.23%. Therefore, the loss of information caused by the constraint is negligible. Now, we consider constraints on the state variables. First, we assume that the changes in C1 (t) in two consecutive discrete time instants of tg are constrained to a region [− C1 , + C1 ]. The design problem includes (12) and the constraints: C1 (ti ) − C1 (ti−1 ) ≤ + C1 ,

ti ∈ t g , i ≥ 1,

(16a)

C1 (ti ) − C1 (ti−1 ) ≥ − C1 ,

ti ∈ t g , i ≥ 1.

(16b)

For demonstrating the treatment of constraints involving state variables let us consider the discretization scheme detailed in

Appendix A. Its application to (16a) produces a set of constraints K  k=0

LK,k (1.0)C1,j,k −

K 

LK,k (1.0)C1,j−1,k ≤ + C1 ,

j ∈ {2, . . ., E},

k=0

(17) which are included in the optimal design problem. Here, C1,j,k is the approximation of C1 at the kth node of the jth finite element. Eq. (16a) is a special case as it involves only approximations of the state variables at the bounds of the finite elements (ti ∈ tg , i ≥ 1), and consequently the discretization only produces additional constraints for C1 at  k = 1.0, where  k is the normalized time variable in finite elements (see Appendix A). This strategy can be general-

98

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

Table 6 DEC: results for 3-COMP batch model with constraints. Constraints (15)

(16)

Optimal design 

constr,u

⎛

=

0.0500,

0.5000,

⎜ [0.0, 0.5[, ⎜ ⎜ ⎜ 0.1536, constr, C1  =⎜ ⎜ [6.0, 6.5[, ⎜ ⎝ 0.4101,

⎛ (18)

0.1948,

0.3948,

0.5000

| 1.0701

| 0.0 0.3739,

0.1949,

0.1061

[0.5, 1.0[,

[1.0, 1.5[,

[1.5, 2.0[,

0.2834,

0.3436,

0.3765

[6.5, 7.0[,

[7.0, 7.5[,

[7.5, 8.0[,

0.4185,

0.5000

| 0.4290

| 0.0 [8.5, 9.0[, [9.0, 9.5[, [9.5, 10.0[ 0.0500, 0.1689, 0.1739, 0.5000



0.0500,

⎟ ⎟ ⎟ 0.3971, ⎟ ⎟ [8.0, 8.5[, ⎟ ⎟ ⎠

CPU time (s)

−9.3032

188.34

−10.2221

207.32

−11.3211

270.92

[2.0, 6.0[,

0.4884,

⎜ [0.0, 1.0[, ⎝

[1.0, 1.5[,

[1.5, 2.0[,

[2.0, 4.0[,

[4.0, 4.5[,

0.0500,

0.4751,

0.3189,

0.2735

| 0.3938

[4.5, 8.5[,

[8.5, 9.0[,

[9.0, 9.5[,

[9.5, 10.0],

| 0.0

 constr,C1 = ⎜

⎞

Optimum

⎞ ⎟ ⎟ ⎠

The first line of  is for u1 (t) and the second for time intervals where it is constant. The last column is for C1in and t0 .

Fig. 9. DEC: dynamic optimal design of experiments for 3-COMP model with constraints (16) on C1 variation: (a) optimal profile of actions; (b) dynamics of the concentrations in the compartments.

ized to constraints holding in the whole time domain which lead to additional constraints at the complete set of nodes of all finite elements, a case addressed next. Consider C1 = 0.05 g/L and call the optimal design  constr, C1 . There are significant differences from the behavior of the unconstrained design in Fig. 7 and the constrained design in Fig. 9. A comparison shows a U-shaped profile of values for u1 (t), rather than one which is initially flat before increasing towards the end of the run, see Fig. 9(a). Consequently, a flatter profile of C1 (t) is obtained as demonstrated in Fig. 9(b). Because of the tight constraint on changes in the value of C1 (t), the magnitude of the changes is small and the values of the concentration in Fig. 9(b), away from a similar minimum value, are less than those for the unconstrained design in Fig. 7(b). The design  constr, C1 is presented in Table 6 (second line) and its efficiency, relative to that of the unconstrained design, is 85.14%. Now, we consider an example to demonstrate the saturation of the state variables; we call saturation to any occurrence where lower (or upper) limits of the variables are reached. For such a purpose we impose upper and lower bounds on C1 (t) over the horizon of the experiment, i.e., the concentration in the first compartment is limited by

C1 (t) ≤ C1max ,

t ∈ [0, H],

(18a)

C1min ,

t ∈ [0, H],

(18b)

C1 (t) ≥

where C1min and C1max are the minimum and maximum values allowed for C1 , respectively. Here, the constraints hold for all the time instants of the discretization grid as well as for the nodes of all the finite elements, and we set C1min to 0.1 g/L and C1max to 0.5 g/L. This design is denoted by  constr,C1 . Fig. 10(a) presents the behavior of the control variable. The increased magnitude of the resulting changes compared to those of the optimal design including the constraints (16) is evident. Fig. 10(b) presents the dynamics of the concentrations in the compartments, and we note that in the interval [1,2] h, C1 (t) is at its minimum which demonstrates the ability of the discretization scheme to efficiently capture the saturation of state variables. The design  constr,C1 appears in Table 6 (third line) and its efficiency, computed as before, is 70.89%. 6. Conclusions Our numerical procedure makes possible the calculation of optimal designs with several kinds of constraints. We find that such a capability is necessary for almost all the problems we have investigated in this work.

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

99

Fig. 10. DEC: dynamic optimal design of experiments for 3-COMP model with constraints (18) on C1 : (a) optimal profile of actions; (b) dynamics of the concentrations in the compartments.

For the static designs of Section 4.1 the only constraints were on the design region, which are standard, and on the maximum concentration of C1 . For the 2-COMP model this latter constraint became binding and the maximum value of u1 (t) was 0.4810 g/(L h), less than the boundary value of 0.5 g/(L h). This constraint did not become active for the 3-COMP model. For the SE&I designs of Section 4.2 one value, of u1 (t) had to be chosen to again ensure that the value of C1 at the end of the experimental run was no more than 1.5 g/L. In the dynamic experiments setup of Section 4.3 the only constraint that is active is that from the design region for the value of u1 (t), which goes from its lowest to its highest value. On the other hand, the three designs of Section 5 involve constraints on the rate of change of either u1 (t) or C1 (t) or in the values of C1 (t). Typically, multi-experiment designs are naturally more informative than single experiment ones. We have given timings of our program in all tables of designs, but have not discussed them. The times to find the optimal design range from around 3–16 min. These times arise because we are interested in finding global optimal designs and must therefore use a global optimizer. Most of the computation time is spent in certifying global optimality. The solver stops when either 100 successive iterations yield no major improvement or 3000 solutions are found. As the timings show, finding designs for 3-COMP models is more challenging, because we have more decision variables and a FIM of higher dimension. Interestingly, a comparison of the timings in Table 5 for DE, the dynamic experiment, and those in Table 6 for the constrained problem DEC, shows that the introduction of constraints leads to appreciable reduction in computational times because the number of degrees of freedom of the later design problem is lower. In our paper we have found D-optimal designs for a specific model with continuous inputs and a response measured over time. The general class of problems to which our method can be applied is much more general than that introduced by Box and Lucas [41] of finding optimum measuring times for the response of batch reactions, a formulation that has been followed in the majority of papers in the statistical literature on optimal design for nonlinear models. Our formulation seems more appropriate for medical and pharmaceutical experiments where treatments and measurements are repeated over time. In particular, the time profile of drugs administered by catheter is similar to that for our dynamic experiments, with changes in level at the discretion of medical staff. Finally, it is noteworthy that compartmental models are a broad class of models with many applications in pharmacokinetics and

pharmacodynamics [30] and our approach allows addressing it successfully. Expectably, it can also be used for identifying more complex models such as those representing the glyceride transesterification kinetics [42] or steam explosion kinetics [43].

Conflict of interest None declared.

Appendix A. Numerical strategy for solving the ODE system We now describe the numerical approach used for solving the dynamic optimization problem described in Section 1. A simultaneous strategy based on orthogonal collocation on finite elements is adopted, where the ODEs are discretized in time, t, and polynomial approximations are used to describe the solution, resulting in a nonlinear system of algebraic equations. Then, the maximization of the logarithm of the determinant of the FIM (D-optimality criterion) is reformulated as a nonlinear programming problem in the compact decision domain, and solved “all-at-once”. In this formulation, first-order derivatives of the variables describing the process dynamics and the sensitivities are expressed as linear combinations of the solutions at the discretization nodes of the finite elements, that correspond to the time instants at which the system is sampled. The duration of the experiments is divided into E unequal finite elements, and a polynomial approximation with order K over a finite element is constructed for Ci (t) as well as for the sensitivities  i,j (t). Here, K is the degree of the approximation of the solution in each finite element, and K + 1 the number of collocation points. To demonstrate the discretization procedure let us consider a single ODE, dy = f (y, t), dt

(A.1a)

y(0) = y0 ,

(A.1b)

where y is the variable to be integrated in a closed interval, f(y, t) is a continuous function and y0 is the initial solution. Then, the approximation of the variable y(t), y˜ (t), is constructed using combinations of orthogonal polynomials at pre-defined local points. In our frame-

100

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101

work y(t) ∈ {Ci (t),  i,j (t)}, i ∈ [ns ], j ∈ [n ]. The discretization allows representing a time instant t located in the finite element j as

⎫ ⎪ ⎪ ⎬

t = tj−1 + hj  y˜ j (t) =

K  k=0

LK,k ()yj,k ⎪ ⎪ ⎭





t ∈ ti,j−1 , ti,j ,

 ∈ [0, 1] ,

(A.2)

[4]

[5] [6]

K

Where LK,k () = l=0,l =/ k ( − l ) / (k − l ) are Kth-order Lagrange interpolation polynomials of the scaled variable . Let the collocation points being  l , l ∈ {0, . . ., K}, where  0 = 0,  l−1 <  l , l ∈ [K], hj is the length of element j, and we use yj,l to compactly designate y(tj,l ), the solution at the lth collocation node of the jth finite element, with tj,l = tj−1,K + hj ×  l . Here, third-order Lagrange interpolation polynomials are used for convenience, i.e., we set K = 3. Substituting the polynomial approximation of y˜ (t) in (A.1), the collocation equations approximating the solution of the ODEs at the collocation points  l , l ∈ [K] can be written: K 

L˙ K,k (l )yj,k = hj × f

 K 



LK,k (l )yj,k , tj,k

l ∈ [K], j ∈ [E].

[7] [8]

[9] [10] [11] [12]

[13]

k=0

k=0

(A.3) Here, L˙ K,k (l ) = dLK,k (l )/d, and f(•) is the right hand function in (A.1a). The interpolation points  l were chosen judiciously so they provide accurate approximations of y˜ (t). For this purpose, the roots of the shifted fourth order Radau polynomials were selected ( 1 = 0.088588,  2 = 0.409467,  3 = 0.787659, and  4 = 1.00000 ≡ ( 0 = 0.00000)) as they assure an exact approximation of y˜ (t) when it is locally described by a general polynomial of order lower than or equal to 2 × K [44,p. 291]. Additionally, constraints are introduced to impose the continuity of y˜ (t) between the finite elements, so that the solution at the first node ( 0 = 0.00000) in element j + 1 coincides with that at the last node ( 4 = 1.00000) of element j; yj+1,0 =

K 

[14]

[15]

[16]

[17] [18]

[19]

[20]

LK,k (1.0)yj,k ,

j ∈ {2, . . ., E}.

(A.4)

[21]

k=0

No constraints are imposed to derivatives d˜y(t)/dt at the bounds of the finite elements since that in result of the sudden “step changes” of the input variables there, the time derivatives may also be discontinuous. Finally, at tUP the approximation of y˜ (t) follows K 

L˙ K,k (1.0)yE,K = hE × f

k=0

 K 

[22]

[23] [24]



LK,k (1.0)yE,k , tE,K

(A.5)

[25]

k=0

and the initial condition reads y1,0 = y0 ,

[26]

(A.6)

where yE,K is the solution of y(t) at tUP , the final integration time, y0 is the initial condition for variable y, and y1,0 is the solution of y(t) at t = 0. The solution of the set of algebraic equations (A.3)–(A.6) is the solution of the original ODE at the discretization nodes. References [1] S. Rosenbaum, Basic Pharmacokinetics and Pharmacodynamics: An Integrated Textbook and Computer Simulations, Wiley, New Jersey, 2012 https://books. google.pt/books?id=HY9iUD2DrrEC. [2] D.R. Brocks, R. Mehvar, Rate and extent of drug accumulation after multiple dosing revisited, Clin. Pharmacokinet. 49 (2010) 421–438, http://dx.doi.org/ 10.2165/11531190-000000000-00000. [3] J.G. Wagner, G.J. Szpunar, J.J. Ferry, Michaelis–Menten elimination kinetics: areas under curves, steady-state concentrations, and clearances for

[27]

[28]

[29] [30] [31]

[32] [33]

compartment models with different types of input, Biopharm. Drug Dispos. 6 (1985) 177–200. A.C. Atkinson, B. Bogacka, Compound and other optimum designs for systems of nonlinear differential equations arising in chemical kinetics, Chemometr. Intell. Lab. Syst. 61 (2002) 17–33, http://dx.doi.org/10.1016/S01697439(01)00173-3. A.C. Atkinson, V.V. Fedorov, The design of experiments for discriminating between two rival models, Biometrika 62 (1975) 57–70. ´ D. Ucinski, B. Bogacka, T-optimum designs for multiresponse dynamic heteroscedastic models, in: A. Di Bucchianico, H. Läuter, H.P. Wynn (Eds.), mODA 7 – Advances in Model-Oriented Design and Analysis, Physica-Verlag, Heidelberg, 2004, pp. 191–199. G. Franceschini, S. Macchietto, Model-based design of experiments: state of the art, Chem. Eng. Sci. 63 (2008) 4846–4872. F. Galvanin, M. Barolo, S. Macchietto, F. Bezzo, Optimal design of clinical tests for the identification of physiological models of type 1 diabetes in the presence of model mismatch, Med. Biol. Eng. Comput. 49 (2011) 263–277, http://dx.doi.org/10.1007/s11517-010-0717-8. ´ D. Ucinski, Optimal Measurement Methods for Distributed Parameter System Identification, CRC Press, Boca Raton, 2005. L. Pronzato, Optimal experimental design and some related control problems, Automatica 44 (2008) 303–325. A.E. Bryson, Y.C. Ho, Applied Optimal Control, Blaisdell, New York, 1969. V.S. Vassiliadis, R.W.H. Sargent, C.C. Pantelides, Solution of a class of multistage dynamic optimization problems. 1. Problems with path constraints, Ind. Eng. Chem. Res. 33 (1994) 2111–2122. M. Diehl, H.G. Bock, J.P. Schl”oder, A real-time iteration scheme for nonlinear optimization in optimal feedback control, SIAM J. Control Optim. 43 (2005) 1714–1736, http://dx.doi.org/10.1137/S0363012902400713. L.T. Biegler, S.L. Campbell, V. Mehrmann, Control and Optimization with Differential-Algebraic Constraints, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2012. V.S. Vassiliadis, R.W.H. Sargent, C.C. Pantelides, Solution of a class of multistage dynamic optimization problems. 2. Problems with path constraints, Ind. Eng. Chem. Res. 33 (1994) 2123–2133. E. Balsa-Canto, D. Henriques, A. Gábor, J.R. Banga, AMIGO2, a toolbox for dynamic modeling, optimization and control in systems biology, Bioinformatics 32 (2016) 3357–3359, http://dx.doi.org/10.1093/ bioinformatics/btw411. I. Bauer, H.G. Bock, S. Körkel, J.P. Schlöder, Numerical methods for optimum experimental design in DAE systems, J. Comput. Appl. Math. 120 (2000) 1–25. S. Körkel, E. Kostina, H.G. Bock, J.P. Schlöder, Numerical methods for optimal control problems in design of robust optimal experiments for nonlinear dynamic processes, Optim. Methods Softw. 19 (2004) 327–338, http://dx.doi. org/10.1080/10556780410001683078. D. Telen, B. Houska, F. Logist, E. Van Derlinden, M. Diehl, J. Van Impe, Optimal experiment design under process noise using Riccati differential equations, J. Process Control 23 (2013) 613–629. G. Carey, B.A. Finlayson, Orthogonal collocation on finite elements, Chem. Eng. Sci. 30 (1975) 587–596, http://dx.doi.org/10.1016/0009-2509(75)80031-5. M.D. Hoang, T. Barz, V.A. Merchan, L.T. Biegler, H. Arellano-Garcia, Simultaneous solution approach to model-based experimental design, AIChE J. 59 (2013) 4169–4183, http://dx.doi.org/10.1002/aic.14145. B.P.M. Duarte, A.C. Atkinson, Optimal designs of experiments for non-isothermal kinetic rates: analysis of different strategies, Optim. Eng. (2018), http://dx.doi.org/10.1007/s11081-018-9415-4 (in press). P. Valko, S. Vajda, An extended ODE solver for sensitivity calculations, Comput. Chem. 8 (1984) 255–271. F. Galvanin, F. Bezzo, Advanced techniques for the optimal design of experiments in pharmacokinetics, in: Quantitative Systems Pharmacology – Models and Model-Based Systems with Applications, Elsevier, 2018, pp. 65–83 (Chapter 3). D. Telen, F. Logist, R. Quirynen, B. Houska, M. Diehl, J. Van Impe, Optimal experiment design for nonlinear dynamic (bio)chemical systems using sequential semidefinite programming, AIChE J. 60 (2014) 1728–1739, http:// dx.doi.org/10.1002/aic.14389. H. Banks, S. Hu, W. Thompson, Modeling and Inverse Problems in the Presence of Uncertainty, Chapman and Hall/CRC, New York, 2014. A. Saltelli, S. Tarantola, K.P.-S. Chan, A quantitative model-independent method for global sensitivity analysis of model output, Technometrics 41 (1999) 39–56, http://dx.doi.org/10.1080/00401706.1999.10485594. T. Maly, L.R. Petzold, Numerical methods and software for sensitivity analysis of differential-algebraic systems, Appl. Numer. Math. 20 (1996) 57–79, http:// dx.doi.org/10.1016/0168-9274(95)00117-4. A.C. Atkinson, A.N. Donev, R.D. Tobias, Optimum Experimental Designs, with SAS, Oxford University Press, Oxford, 2007. V.V. Fedorov, S.L. Leonov, Optimal Design for Nonlinear Response Models, Chapman and Hall/CRC Press, Boca Raton, 2014. A. Atkinson, V. Fedorov, A. Herzberg, R. Zhang, Elemental information matrices and optimal experimental design for generalized regression models, J. Stat. Plann. Inference 144 (2014) 81–91, http://dx.doi.org/10.1016/j. jspi2012.09.012. V.V. Fedorov, The design of experiments in the multiresponse case, Theory Probabil. Appl. 16 (1971) 323–332. J.J. Dongarra, C.B. Moler, J.R. Bunch, G.W. Stewart, LINPACK User’s Guide, 1979, Philadelphia, PA.

B.P.M. Duarte et al. / Journal of Process Control 83 (2019) 88–101 [34] G. Golub, C. Van Loan, Matrix Computations, Johns Hopkins Studies in the Mathematical Sciences, Johns Hopkins University Press, 2013 https://books. google.pt/books?id=X5YfsuCWpxMC. [35] GAMS Development Corporation, GAMS – A User’s Guide, GAMS Release 24.2.1, 2013, Washington, DC, USA. [36] Z. Ugray, L. Lasdon, J. Plummer, F. Glover, J. Kelly, R. Martí, A multistart scatter search heuristic for smooth NLP and MINLP problems, in: Metaheuristic Optimization via Memory and Evolution, Springer, 2005, pp. 25–51. [37] A. Drud, CONOPT: a GRG code for large sparse dynamic nonlinear optimization problems, Math. Programm. 31 (1985) 153–191. [38] W. Chen, K. Wang, Z. Shao, L.T. Biegler, Direct transcription with moving finite elements, in: Control and Optimization with Differential-Algebraic Constraints, SIAM, 2012, pp. 233–252, http://dx.doi.org/10.1137/ 9781611972252.ch11 (Chapter 11). [39] A. Akintola, S. Jansen, R. Wilde, G. Hultzer, R. Rodenburg, D. van Heemst, A simple and versatile method for frequent 24 h blood sample collection in

[40]

[41] [42]

[43]

[44]

101

healthy older adults, MethodsX 2 (2015) 33–38, http://dx.doi.org/10.1016/j. mex.2014.12.003. M.B. Zarrop, Optimal Experimental Design for Dynamic System Identification: Lecture Notes in Control and Information Sciences 21, Springer-Verlag, New York, 1979. G.E. Box, H. Lucas, Design of experiments in non-linear situations, Biometrika 46 (1959) 77–90. J.F.O. Granjo, B.P.M. Duarte, N.M.C. Oliveira, Systematic development of kinetic models for the glyceride transesterification reaction via alkaline catalysis, Ind. Eng. Chem. Res. 57 (2018) 9903–9914, http://dx.doi.org/10. 1021/acs.iecr.7b05328. R. Vargas, A. Vecchietti, Modeling the thermochemical pretreatment of eucalyptus globulus for bioethanol production, Ind. Eng. Chem. Res. 57 (2018) 12458–12467, http://dx.doi.org/10.1021/acs.iecr.8b02706. L.T. Biegler, Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2010.