Robust Optimization for A Nonlinear Switched Time-Delay System with Noisy Output Measurements Using Hybrid Optimization Algorithm

Robust Optimization for A Nonlinear Switched Time-Delay System with Noisy Output Measurements Using Hybrid Optimization Algorithm

Accepted Manuscript

Robust Optimization for A Nonlinear Switched Time-Delay System with Noisy Output Measurements Using Hybrid Optimization Algorithm Jinlong Yuan, Jun Xie, Chongyang Liu, Kok Lay Teo, Ming Huang, Houming Fan, Enmin Feng, Zhilong Xiu PII: DOI: Reference:

S0016-0032(19)30484-3 https://doi.org/10.1016/j.jfranklin.2019.06.037 FI 4026

To appear in:

Journal of the Franklin Institute

Received date: Revised date: Accepted date:

19 September 2018 1 March 2019 23 June 2019

Please cite this article as: Jinlong Yuan, Jun Xie, Chongyang Liu, Kok Lay Teo, Ming Huang, Houming Fan, Enmin Feng, Zhilong Xiu, Robust Optimization for A Nonlinear Switched Time-Delay System with Noisy Output Measurements Using Hybrid Optimization Algorithm, Journal of the Franklin Institute (2019), doi: https://doi.org/10.1016/j.jfranklin.2019.06.037

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

00 (2019) 1–??

CR IP T

Robust Optimization for A Nonlinear Switched Time-Delay System with Noisy Output Measurements Using Hybrid Optimization Algorithm Jinlong Yuana,b , Jun Xiec , Chongyang Liud , Kok Lay Teoe,f , Ming Huanga , Houming Fanb,∗, Enmin Fengg , Zhilong Xiuh a Department

of Mathematics, School of Science, Dalian Maritime University, Dalian 116026, Liaoning, P.R. China of Transportation Engineering, Dalian Maritime University, Dalian 116026, Liaoning, P.R. China and Research Office of Mathematics, Department of Basics, PLA Dalian Naval Academy, Dalian 116018, Liaoning, P.R. China d School of Mathematics and Information Science, Shandong Institute of Business and Technology, Yantai 264005, PR China e Department of Mathematics and Statistics, Curtin University, Perth, WA 6102, Australia f The Coordinated Innovation Center for Computable Modeling in Management Science, Tianjin University of Finance and Economics, Tianjin 300222, China g School of Mathematical Sciences, Dalian University of Technology, Dalian 116024, Liaoning, P.R. China h School of Life Science and Biotechnology, Dalian University of Technology, Dalian 116024, Liaoning, P.R. China

M

AN US

b College

c Teaching

Abstract

CE

PT

ED

This paper considers the problem of using noisy output data to estimate unknown switching times and unknown system parameters arising in 1,3-PD batch production process of glycerol induced by Klebsiella pneumonia. We formulate the problem as a robust optimization problem in which the unknown quantities are decision variables to be chosen optimally, with the cost function penalizing the mean and variance of the relative error between the output of the system and the measured actual noisy system output. This problem is governed by a switched time-delay system subject to continuous state inequality constraints arising from engineering specifications. A new time-scaling transformation and constraint transcription technique are used then to convert this resulting problem into a sequence of approximate subproblems. A hybrid optimization algorithm combined genetic algorithm with gradient-based method is developed to solve these resulting subproblems. Finally, we explore the correctness of the optimal switching times and system parameters obtained, as well as the effectiveness of the proposed algorithm. c 2019 Published by Elsevier Ltd.

AC

Keywords: Nonlinear programming, Robust optimization, Switched time-delay system, Noisy output data, Hybrid optimization algorithm

1. Introduction

1,3-Propanediol (1,3-PD), a valuable bifunctional molecule, is a monomer widely used in the production of polymers, such as polyesters and polyurethanes [1]. Generally, 1,3-PD is produced by chemical synthesis or the bioconversion of glycerol induced by Klebsiella pneumonia (K. pneumonia) [2]. The latter approach is now becoming increasingly attractive in industry because it is environmentally friendly and uses a cheap renewable feedstock [3]. At ∗ Corresponding

author Email addresses: [email protected]; [email protected] (Jinlong Yuan), [email protected] (Houming Fan)

1

ACCEPTED MANUSCRIPT

2

/ 00 (2019) 1–??

AC

CE

PT

ED

M

AN US

CR IP T

present, the fermentation process can be divided into three categories: batch, continuous and fed-batch fermentations [4, 5]. In batch fermentation, bacteria and substrate are presented at the beginning of the process, and nothing is added or removed from the fermentor during the process. In continuous fermentation, fresh medium is added continuously to replenish consumed substrate while old medium is removed during the reaction. Fed-batch fermentation is a mixture of the batch and continuous cultures. This paper focuses on batch fermentation because it lays a foundation for continuous and fed-batch fermentations (i.e. it is a part of fed-batch or continuous fermentations) [6]. Consequently, batch fermentation has seen a recent vigorous surge of activity in several related research areas over the last two decades. These include robust identification [7], robust optimal control problems [8], robust bi-objective optimal control [9], sensitivity analysis [10], stochastic system [11], strong stability [12], two-stage system [13] and switched system [14]. However, in aforementioned studies, the presence of time-delay is ignored. In real world, many practical problems naturally contain time delays and they must not be ignored; see, for example, [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] and references therein. During the fermentation process, biomass and glycerol are added to the reactor. However, the reaction will not take place instantaneously. It needs to go through some processes: (i) biomass and glycerol need to be mixed thoroughly; and (ii) glycerol will move through cell membranes from outside to inside. These processes do take certain time to accomplish. Thus, time-delay system is an ideal tool to characterize the processes mentioned above [26]. Recently, Liu et al. [27] and Yuan et al. [28] studied parameter identification for a nonlinear time-delay system in microbial batch fermentation. Experimental design suboptimization for the enzyme-catalytic nonlinear time-delay system in microbial batch culture is studied in [29]. Furthermore, robust optimization for nonlinear time-delay dynamical system of dha regulon with cost sensitivity constraint in batch culture is investigated in [30]. For all the papers mentioned above, a major drawback is that they do not take into account the presence of noise in the output data. That is, it is assumed that the real system output data (measured through experiments) are obtained exactly. Clearly, this is an idealistic assumption. In practice, identical experiments are unlikely to yield exactly the same results. There will always be some randomness due to some disturbances or measurement errors [31]. For noisy output measurements in batch fermentation, it is in its infancy. This motivates the work in this paper. In this paper, a nonlinear switched time-delay system involving unknown switching times and unknown system parameters is proposed to describe the batch fermentation of glycerol induced by K. pneumonia to 1,3-PD. Our goal is to estimate the unknown quantities by using noisy output data. To this end, we present a robust optimization problem governed by the switched time-delay system, where the cost function is a weighted sum of the expectation and variance of the relative error between the output of the switched time-delay system and the measured noisy actual system output. The unknown quantities are regarded as decision variables to be chosen optimally, subject to continuous state inequality constraints. Then, this problem is converted into a sequence of approximate subproblems by using a novel combination of an new time-scaling transformation [32], constraint transcription technique [33]. Furthermore, an algorithm combined genetic algorithm and gradient-based optimization is developed to solve each of these approximate subproblems. Finally, numerical results show that the nonlinear switched time-delay system with the optimal switching times and system parameters can reasonably describe the process of batch fermentation. The remainder of this paper is organized as follows. In Section 2, we give the nonlinear switched time-delay. In Section 3, we present a robust optimization problem. In Section 4, numerical solution methods are developed to solve this optimization problem. In Section 5, numerical results are presented. In Section 6, we draw some useful conclusions and indicate the direction for future research. 2. Nonlinear switched time-delay system According to the actual experiment process, we make the following assumptions: A1: No medium is pumped inside and outside the bioreactor in the process of batch fermentation. A2: The concentrations of reactants are uniformly distributed in the reactor. A3: There are no state jumps at the switching times. Nomenclature n o • In denotes the set 1, 2, . . . , n . • R+ denotes the set of nonnegative real numbers.

2

ACCEPTED MANUSCRIPT

3

/ 00 (2019) 1–??

• R− denotes the set of negative real numbers. • R denotes the set of real numbers.

• AT denotes the transposition of matrix or vector A.

• C 1 (R− , R5+ ) denotes the Banach space of continuously differentiable functions mapping from R− into R5+ .

• h denotes a given time-delay, which is given in [27], for biomass in the model of the fermentation process.

• [0, t f ] ⊂ R+ denotes the interval of reaction time. • N denotes the total number of the experiments.

CR IP T

• t f denotes a given terminal time for the fermentation process.

• Ωl , l ∈ IN , denote the set of given characteristic time points (two or more discrete time points) during the culture process of the lth experiment on the time horizon of t ∈ [0, t f ]. e l , l ∈ IN , denote the set of given characteristic time points (two or more discrete time points) during the culture • Ω process of the lth experiment on the time horizon of s ∈ [0, 3]. e l |, l ∈ IN , and |Ωl |, l ∈ IN , denote the numbers of elements of the sets Ωl and Ω e l , respectively. • |Ω

AN US

l l T • x0l := (x01 , . . . , x05 ) ∈ R5 , l ∈ IN , denote the initial state vectors of the lth experiment.

• xi,l (t) = (x1i,l (t), . . . , x5i,l (t))T ∈ R5 , l ∈ IN , i ∈ I3 , denote the state vectors whose components are, respectively, the concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, acetate and ethanol at time t ∈ [τi−1 , τi ] of the lth experiment. • xi−1,l (τi−1 ), l ∈ IN , i = 2, 3, denotes both the end state vectors of the (i − 1)th switched stage and the initial state vectors of the ith switched stage of the lth experiment. • µim , i ∈ I3 , are the maximum specific growth rates of the ith switched stage.

M

• k1i , i ∈ I3 , are the Monod saturation constants of the ith switched stage.

• mi2 , i ∈ I3 , are the maintenance terms of substrate consumption under substrate-limited conditions of the ith switched stage.

ED

• Y2i , i ∈ I3 , are the maximum growth yields of the ith switched stage.

• mij , j = 3, 4, 5, i ∈ I3 , are the maintenance terms of 1,3-PD, acetate and ethanol under substrate-limited conditions of the ith switched stage.

PT

• Y ij , j = 3, 4, 5, i ∈ I3 , are the maximum 1,3-PD, acetate and ethanol yields of the ith switched stage.

CE

In the batch fermentation, there are three phases, i.e., lag phase, exponential growth phase and stationary phase [41]. During the lag phase, the cells will take time to adjust and the exponential phase is characterized by cell doubling. During the stationary phase, growth will slow down and eventually stop as nutrients become depleted or inhibitory metabolites build up in this phase [41]. With this in mind, there are two different switching times τ1 and τ2 such that 0 = τ0 < a < τ1 < b < c < τ2 < d < τ3 = t f .

(1)

AC

Let τ := (τ1 , τ2 )T and let D denote the set of all τ satisfied condition (1). The system parameter vectors are also to be chosen optimally. They are defined as: ζ i = (µim , k1i , mi2 , Y2i , mi3 , Y3i , mi4 , Y4i , mi5 , Y5i )T ∈ R10 , i ∈ I3 .

Define

Z :=

10 Y [z j∗ , z j∗ ],

(2)

j=1

where z j∗ = z0j − 0.5|z0j |, z j∗ = z0j + 0.5|z0j |, j ∈ I10 , are, respectively, the lower and upper bounds of kinetic parameters, and z0 ∈ R10 + is the nominal value of the system parameters. These values are obtained from a previously performed 3

ACCEPTED MANUSCRIPT

4

/ 00 (2019) 1–??

where φl (t) ∈ C 1 (R− , R5+ ), l ∈ IN , are given initial functions; and

CR IP T

parameter estimation procedure in [27]. Let Z3 be the set of all those ζ := ((ζ 1 )T , (ζ 2 )T , (ζ 3 )T )T ∈ R30 such that ζ i ∈ Z, i ∈ I3 . Accordingly, any pair (τ, ζ) ∈ D × Z3 is called an admissible switching-parameter pair. For batch fermentation, the mass balance relationships for biomass, substrate and reaction products can be expressed by the following nonlinear switched time-delay system:  i,l  x˙ (t) = f (i,l (xi,l (t), xi,l (t − h), ζ i ), t ∈ [τi−1 , τi ], i ∈ I3 , l ∈ IN ,      x0l , t = 0, i = 1, l ∈ IN ,  i,l x (t) = (3)  i−1,l   x (t), t = τi−1 , i = 2, 3, l ∈ IN ,     x1,l (t) = φl (t), t ≤ 0, l ∈ IN , f i,l (xi,l (t), xi,l (t − h), ζ i ) = ( f1i,l (xi,l (t), xi,l (t − h), ζ i ), . . . , f5i,l (xi,l (t), xi,l (t − h), ζ i ))T , i ∈ I3 , l ∈ IN ,  i,l i,l  f (x (t), xi,l (t − h), ζ i ) = µi,l (t)x1i,l (t − h),     1i,l i,l i,l f2 (x (t), xi,l (t − h), ζ i ) = −qi,l  2 (t)x1 (t − h),    i,l i,l i,l  f (xi,l (t), xi,l (t − h), ζ i ) = q (t)x (t − h), j∈ {3, 4, 5}. j j 1

(4)

AN US

In particular,

M

The specific substrate consumption rate and the specific product formation rates can be expressed as follows:  5  xi,lj (t)  x2i,l (t) Y   i i,l  µ (t) = µ (1 − ),  m i,l   x∗j  x2 (t) + k1i j=2       i,l µi,l (t) (5)  q2 (t) = mi2 + ,    Y2i       qi,lj (t) = −mij + µi,l (t)Y ij , j∈ {3, 4},      qi,l (t) = mi + µi,l (t)Y i . 5 5 5

ED

It should be noted that there exist critical concentrations, outside which cells cease to grow, of biomass, glycerol, 1,3-PD, acetate and ethanol. As a result, it is biologically meaningful to restrict the concentrations of biomass, glycerol and products to be in a set W defined by

PT

W := [x∗ , x∗ ] =

5 Y [x j∗ , x∗j ] ⊂ R5+ ,

(6)

j=1

CE

where x∗ = [0.0001, 0.1, 0, 0, 0]T , x∗ = [15, 2039, 939.5, 1026, 360.9]T (see [27]). 3. Robust optimization problem

AC

Let xl (·|τ, ζ) denote the solution of system (3) corresponding to each switching-parameter pair (τ, ζ) ∈ D × Z3 . That is, xl (·|τ, ζ) satisfies Z t    1,l l 1 l  x (t; x , τ, ζ ) = x + f 1,l (x1,l (λ), x1,l (λ − h), ζ 1 )dλ, t ∈ [0, τ1 ],   0 0   0  Z t     2,l 1,l 2 1,l l x (t; x (τ ), τ, ζ ) = x (τ ) + f 2,l (x2,l (λ), x2,l (λ − h), ζ 2 )dλ, t ∈ (τ1 , τ2 ], x (t|τ, ζ)=  (7) 1 1     Zτ1t     3,l 2,l 3 2,l   f 3,l (x3,l (λ), x3,l (λ − h), ζ 3 )dλ, t ∈ (τ2 , τ3 ],  x (t; x (τ2 ), τ, ζ ) = x (τ2 ) + τ2

and xl (t|τ, ζ) = φl (t), ∀t ≤ 0, l ∈ IN . Let yl (tk ) := (yl1 (tk ), yl2 (tk ), yl3 (tk ))T ∈ R3 denote the lth experiment data of the concentration of the first three extracellular substances at the time points tk ∈ Ωl , k ∈ I|Ωl | , l ∈ IN , with initial state x0l and φl (t), ∀t ≤ 0, l ∈ IN . 4

ACCEPTED MANUSCRIPT

5

/ 00 (2019) 1–??

In [7], the output measurements yl (tk ), tk ∈ Ωl , k ∈ I|Ωl | , l ∈ IN , are assumed to be obtained exactly. Clearly, this assumption is unrealistic. In reality, the output measurements are always inexact, caused by factors, such as measurement errors. Therefore, in this paper, we assume that yl (tk ), tk ∈ Ωl , k ∈ I|Ωl | , l ∈ IN , are multivariate normally distributed random vectors, where the mean vector (of dimension 3|Ωl |) and the covariance matrix (of dimension 3|Ωl | × 3|Ωl |) are known. To measure the estimation accuracy, we have the following definition: Definition 1. For a switching-parameter pair (τ, ζ) ∈ D × Z3 , the relative error between xlj (tk |τ, ζ), j ∈ I3 , tk ∈ Ωl , k ∈ I|Ωl | , l ∈ IN , and y j (tk ), j ∈ I3 , tk ∈ Ωl , k ∈ I|Ωl | , l ∈ IN , is defined by 5N

1 PN

l=1

|Ωl |

N X 5 X n xl (t |τ, ζ) − yl (t ) o X j k j k 2 l=1 j=1 tk ∈Ωl

x∗j

CR IP T

Re (τ, ζ) :=

.

AN US

Note that Re (τ, ζ) is a function containing random vectors yl (tk ), tk ∈ Ωl , k ∈ I|Ωl | , l ∈ IN . One of our aims is to choose a switching-parameter pair (τ, ζ) ∈ D × Z3 so that the expected value of Re (τ, ζ) is minimized (i.e., the measured system output should be close to the actual system output ”on average”). Also, the variance of Re (τ, ζ) should be minimized to ensure that the optimal estimates for switching-parameter pair (τ, ζ) ∈ D × Z3 , is robust with respect to uncertainties in the output data. Thus, the cost function is defined by n o n o J(τ, ζ) := αE Re (τ, ζ) + (1 − α)D Re (τ, ζ) , (8) n o n o where (τ, ζ) ∈ D × Z3 , E · is the expectation, D · is the variance, and α ∈ [0, 1] is a given weight.

ED

M

Remark 1. The first term of (8) is only regarded as the cost function provided that the distribution of the measurement data is known exactly. Clearly, this is an idealistic assumption. There exists uncertainty in the distribution of measurement data. Thus, it is significant to find the optimal estimates for switching-parameter pair (τ, ζ) ∈ D × Z3 which is robust against the uncertainty. With this in mind, we take the variance of Re (τ, ζ), which measures how far the random error functions are spread out from their expectation value, as the second part of the cost function. The constraints xi,lj (t; xi−1,l (τi−1 ), τ, ζ i ) ∈ W, j ∈ I8 , i ∈ I3 , l ∈ IN , ∀t ∈ [τi−1 , τi ], can be expressed explicitly as follows: gi,lj (xi,l (t; xi−1,l (τi−1 ), τ, ζ i )) := xi,lj (t; xi−1,l (τi−1 ), τ, ζ i ) − x∗j ≤ 0, ∀t ∈ [τi−1 , τi ], j ∈ I5 , i ∈ I3 , l ∈ IN ,

gi,lj+5 (xi,l (t; xi−1,l (τi−1 ), τ, ζ i ))

:= x j∗ −

xi,lj (t; xi−1,l (τi−1 ), τ, ζ i )

≤ 0, ∀t ∈ [τi−1 , τi ], j ∈ I5 , i ∈ I3 , l ∈ IN .

(9a) (9b)

CE

PT

These constraints are also called the path constraints or continuous state inequality constraints. Each of these constraints constitutes an infinite number of restrictions to be satisfied at every time point in the time horizon [33, 35]. Clearly, gi,lj : R → R, i ∈ I3 , l ∈ IN , j ∈ I10 , are continuously differentiable functions. The robust optimization problem governed by system (3) and subject to continuous state inequality constraints can be stated formally as follows:

AC

Problem OP: Given system (3), choose a switching-parameter pair (τ, ζ) ∈ D × Z3 such that J(τ, ζ) is minimized over D × Z3 subject to continuous state inequality constraints (9a) and (9b).

Remark 2. The goal in Problem OP is to minimize both the average relative error and the relative error variance. The relative importance between these two objectives is controlled by the weight α. If α approaches one, then the average relative error is preferentially minimized; if α approaches zero, then the relative error variance is preferentially minimized. In the case of exactly knowing the output distribution, the best option for minimizing the expected relative error is α = 1. However, the measured output distribution are only approximations of the true one due to errors and/or uncertainties in the output distribution, it is indispensable to select α < 1 to guarantee solution robustness.

Theorem 1. Problem OP admits an optimal solution. 5

ACCEPTED MANUSCRIPT

6

/ 00 (2019) 1–??

Proof. From (8), (9a) and (9b), we see that the cost function J(τ, ζ) is continuous in (τ, ζ) ∈ D × Z3 and the constraint functions gi,lj , j ∈ I5 , i ∈ I3 , l ∈ IN and gi,lj+5 , j ∈ I5 , i ∈ I3 , l ∈ IN are continuous in (t, τ, ζ) ∈ [0, t f ] × D × Z3 . Since the set D × Z3 is compact, the proof follows readily.  To proceed, we shall show how to compute the cost function J(τ, ζ) in the following theorem. Theorem 2. For each (τ, ζ) ∈ D × Z3 , J(τ, ζ) =

5N

α PN

N X 5 X ( xl (t |τ, ζ)  X 2 j k

|Ωl |

l=1

x∗j

l=1 j=1 tk ∈Ωl

n ylj (tk ) o xlj (tk |τ, ζ) n ylj (tk ) 2 o) − 2E + +E x∗j x∗j x∗j

1

CR IP T

( X N X N X X X 5 xl1 (t |τ, ζ) 5 X n ylj11 (tk1 ) ylj22 (tk2 ) o 1−α j1 k 1 · Cov , · 4 P l 2 x∗j1 x∗j1 x∗j2 25N 2 ( N l=1 Ω ) l1 =1 l2 =1 tk ∈Ωl tk ∈Ωl j1 =1 j2 =1 2

N XX 5  yl (t )  o) n ylj22 (tk2 ) oi nX j k 2 , − 2E +D x∗j2 x∗j l=1 t ∈Ω j=1

h xlj22 (tk2 |τ, ζ) x∗j2

k

n o where Cov ·, · denotes the covariance. 1 PN

5N

=

l=1

1 PN

5N

l=1

Hence,

l=1 j=1 tk ∈Ωl

5N

n o D Re (τ, ζ)

=

25N 2 (

CE

x∗j

l=1

−2

ylj (tk ) xlj (tk |τ, ζ) x∗j

N X 5 X ( xl (t |τ, ζ)  X 2 j k

|Ωl |

x∗j

l=1 j=1 tk ∈Ωl

n ylj (tk ) 2 o) +E , x∗j

1 PN

−4Cov

AC

N X 5 X ( xl (t |τ, ζ)  X 2 j k

1 PN

PT

and

=

|Ωl |

x∗j

l=1 j=1 tk ∈Ωl

ED

n o E Re (τ, ζ)

|Ωl |

5 X n xl (t |τ, ζ) − yl (t ) o N X X j k j k 2

M

Re (τ, ζ) =

AN US

Proof. The relative error Re (τ, ζ) can be written as follows:

l=1

Ωl )2

(

4D

x∗j

+

 ylj (tk ) 2 ) x∗j

tk1 ∈Ω

.

n ylj (tk ) o xlj (tk |τ, ζ) − 2E x∗j x∗j (11)

N X X 5 yl (t ) xl (t |τ, ζ) o nX j k j k l=1 tk ∈Ωl

j=1

x∗j

x∗j

+D

N X X 5  yl (t )  o nX j k 2 l=1 tk ∈Ωl

j=1

5  yl2 (t )  o) N X X 5 yl1 (t ) xl1 (t |τ, ζ)  X N X X nX j2 k2 2 j1 k1 j1 k1 , . ∗ ∗ x x x∗j2 l l j1 j1 j =1 j =1 l =1 l =1 1

(10)

1

2

tk2 ∈Ω

x∗j

(12)

2

Note that

D

N X X 5 yl (t ) xl (t |τ, ζ) o nX j k j k l=1 tk ∈Ωl

=

j=1

x∗j

N X N X X X 5 X 5 X

l1 =1 l2 =1 tk1 ∈Ωl tk2 ∈Ωl j1 =1 j2 =1

=

x∗j

Cov

n ylj11 (tk1 ) xlj11 (tk1 |τ, ζ) ylj22 (tk2 ) xlj22 (tk2 |τ, ζ) o , x∗j1 x∗j1 x∗j2 x∗j2

N X N X X X 5 X 5 xl1 (t |τ, ζ) X j1 k1 l1 =1 l2 =1 tk1 ∈Ωl tk2 ∈Ωl j1 =1 j2 =1

x∗j1

6

Cov

n ylj11 (tk1 ) ylj22 (tk2 ) o xlj22 (tk2 |τ, ζ) , . x∗j1 x∗j2 x∗j2

(13)

ACCEPTED MANUSCRIPT

7

/ 00 (2019) 1–??

Furthermore, N X X N X X 5 yl1 (t ) xl1 (t |τ, ζ)  X 5  yl2 (t )  o nX j1 k1 j2 k2 2 j1 k 1 Cov , ∗ ∗ x x x∗j2 l l j1 j1 j =1 j =1 l =1 l =1 1

1

tk2 ∈Ω

2

5 X

Substituting (13) and (14) into (12) gives =

1

From the definition of covariance, it follows that

ylj11 (tk1 ) x∗j1

and

n ylj11 (tk1 )  ylj22 (tk2 ) 2 o n ylj11 (tk1 )  ylj22 (tk2 ) 2 o n ylj11 (tk1 ) o n ylj22 (tk2 ) 2 o , = E − E E . x∗j1 x∗j2 x∗j1 x∗j2 x∗j1 x∗j2

ylj22 (tk2 ) x∗j2

=

n ylj1 (tk1 ) o E 1∗ x j1

ED

n ylj1 (tk1 ) ylj2 (tk2 ) o E 1 ∗ 2 ∗ x j1 x j2

Cov

+

PT

n yl2 (tk ) o where E · j2x∗ 2 denotes conditional expectation given j2

AC

CE

n ylj11 (tk1 )  ylj22 (tk2 ) 2 o E x∗j1 x∗j2

n ylj1 (tk1 ) 1

(16)

D

l

n

y j2 (tk2 ) 2

x∗j

y j2 (tk2 ) o l

,

x∗j 1 l y j2 (tk2 ) 2 x∗j 2

2

x∗j

2

o

h ylj22 (tk2 ) x∗j2

n ylj2 (tk2 ) oi −E 2∗ , x j2

. Thus, by the law of total expectation,

2

( n yl1 (t )  yl2 (t )  yl2 (t ) o) j k1 j2 k2 2 j2 k2 = E E 1∗ ∗ x j1 x∗j2 x j2 ( yl2 (t )  n yl1 (t ) yl2 (t ) o) j2 k2 2 j k1 j k2 = E E 1 ∗ 2 ∗ ∗ x j2 x j1 x j2

n ylj11 (tk1 ) o n ylj22 (tk2 ) 2 o = E E + x∗j1 x∗j2 n yl1 (tk ) yl2 (tk ) o Cov j1x∗ 1 , j2x∗ 2 h n yl2 (tk ) 3 o n ylj22 (tk2 ) n ylj22 (tk2 ) 2 ooi j1 j2 j2 2 E − E E . n ylj2 (tk2 ) o x∗j2 x∗j2 x∗j2 2 D x∗ j2

Using the identities

(15)

are joint normal random variables, it follows from Theorem 3.3.1 in [34] that

M

Since

2

n ylj11 (tk1 ) ylj22 (tk2 ) o xlj22 (tk2 |τ, ζ) n ylj11 (tk1 )  ylj22 (tk2 ) 2 oi h Cov , − Cov , x∗j1 x∗j2 x∗j2 x∗j1 x∗j2 N X X 5  yl (t )  o) nX j k 2 . +D x∗j l j=1 l=1 tk ∈Ω

Cov

(14)

( X N X N X X X 5 X 5 xl1 (t |τ, ζ) j1 k1 4 · 2 l 2 x∗j1 25N ( l=1 Ω ) l1 =1 l2 =1 tk ∈Ωl tk ∈Ωl j1 =1 j2 =1 1 PN

AN US

n o D Re (τ, ζ)

2

X X n ylj11 (tk1 )  ylj22 (tk2 ) 2 o xlj11 (tk1 |τ, ζ) Cov , . x∗j1 x∗j1 x∗j2 l1 =1 l2 =1 tk1 ∈Ωl tk2 ∈Ωl j1 =1 j2 =1 5 X

CR IP T

=

tk1 ∈Ω

N X N X

n ylj22 (tk2 ) 2 o n ylj22 (tk2 ) o n ylj22 (tk2 ) o 2 E = D + E x∗j2 x∗j2 x∗j2 7

(17)

ACCEPTED MANUSCRIPT

8

/ 00 (2019) 1–??

and

n ylj22 (tk2 ) 3 o n ylj22 (tk2 ) o n ylj22 (tk2 ) o n ylj22 (tk2 ) o 3 E = 3E D + E , x∗j2 x∗j2 x∗j2 x∗j2

(17) is reduced to

n ylj11 (tk1 )  ylj22 (tk2 ) 2 o n ylj11 (tk1 ) o n ylj22 (tk2 ) 2 o n ylj22 (tk2 ) o n ylj11 (tk1 ) ylj22 (tk2 ) o E =E E − 2E Cov , . x∗j1 x∗j2 x∗j1 x∗j2 x∗j2 x∗j1 x∗j2

Cov

CR IP T

Substituting (18) into (16) yields

n ylj22 (tk2 ) o n ylj11 (tk1 ) ylj22 (tk2 ) o n ylj11 (tk1 )  ylj22 (tk2 ) 2 o , = 2E Cov , . x∗j1 x∗j2 x∗j2 x∗j1 x∗j2

(18)

(19)

From (19), (11), (15) and (8), the proof follows readily.  N X X 5  yl (t )  o nX n ylj (tk ) 2 o j k 2 and D on the right-hand side of equation (10) are indeSince the terms of E ∗ xj x∗j l j=1 l=1 tk ∈Ω

Problem EOP: Given system (3), such that H(τ, ζ)

=

3N

α PN

l=1

|Ωl |

AN US

pendent of τ and ζ, Problem OP is equivalent to the following optimization problem.

N X 3 X ( xl (t |τ, ζ)  X 2 j k l=1 j=1 tk ∈Ωl

x∗j

n ylj (tk ) o xlj (tk |τ, ζ) ) − 2E + x∗j x∗j

( X N X N X X X 3 X 3 xl1 (t |τ, ζ) n ylj11 (tk1 ) ylj22 (tk2 ) o 1−α j1 k1 · Cov , 4 P l 2 x∗j1 x∗j1 x∗j2 9N 2 ( N l=1 Ω ) l1 =1 l2 =1 tk ∈Ωl tk ∈Ωl j1 =1 j2 =1 2

M

1

(20)

ED

n ylj22 (tk2 ) oi) h xlj22 (tk2 |τ, ζ) · − 2E , x∗j2 x∗j2

is minimized over D × Z3 subject to continuous state inequality constraints (9a) and (9b).

PT

4. Numerical solution methods

AC

CE

4.1. A new time-scaling transformation In Problem EOP, the switching times τ1 and τ2 are decision variables to be chosen optimally. It is cumbersome to integrate the state and variational or costate systems numerically when the switching times τ1 and τ2 are decision variables. In [37], the time-scaling transformation is introduced to deal with optimal control problems with variable switching times points. It is now one of the most popular tools to optimize switching times in dynamic optimization problems. This transformation, however, is only applicable to systems without time-delay. In [32, 36], a new timescaling transformation strategy is proposed to optimize switching times in dynamic optimization subproblems. Let ϑi denote the duration of the ith subinterval in (τi−1 , τi ]. That is, ϑi := τi − τi−1 , i ∈ I3 .

(21)

The time-scaling transformation works by introducing a new time variable s ∈ [−h, 3] and relating s to t through the equation   tf , if s = 3,     bsc  X   t = µ(s|ϑ) ¯ = (22) ϑ` + ϑbsc+1 (s − bsc), if s ∈ [i − 1, i), i ∈ I3 ,     `=1    s, if s < 0, 8

ACCEPTED MANUSCRIPT

9

/ 00 (2019) 1–??

where µ¯ is the so-called time-scaling function and b·c denotes the floor function. It is easy to verify that the time-scaling function is non-decreasing, continuous and piecewise linear. Moreover, dµ(s|ϑ) ¯ = ϑbsc+1 , s ∈ (i − 1, i), i ∈ I3 , ds

(23)

where the new decision parameters ϑi , i ∈ I3 , satisfy

i=1

(24)

CR IP T

0 < ϑi ≤ t f , i ∈ I3 , 3 X ϑi = t f .

(25)

Let ϑ := [ϑ1 , ϑ2 , ϑ3 ]T ∈ R3 . Let V denote the set of all duration vectors such that the constraints (24) and (25) are satisfied. Now, since i − (i − 1) = 1 for i ∈ I3 , µ(i|ϑ) ¯ =

i X

ϑ` = τi .

`=1

AN US

Therefore, the time-scaling function maps the fixed integer s = i to the switching time t = τi . Let x˜i,l (s) := xi,l (µ(s|ϑ)), ¯ i ∈ I3 , l ∈ IN . Then, applying the time-scaling transformation t = µ(s|ϑ) ¯ gives the following new equivalent dynamic system:  i,l 

¯  bsc+1 dx (µ(s|ϑ))  ˙˜i,l (s) = d x˜i,l (s) = d xi,l (µ(s|ϑ))  x ¯ = ϑ , s ∈ (i − 1, i], i ∈ I3 , l ∈ IN ,    ds ds dt ( l   x0 , s = 0, i = 1, l ∈ IN , (26)   x˜i,l (s) =  i−1,l   x ˜ (s), s = i, i = 2, 3, l ∈ I , N     x˜1,l (s) = φl (s), s ≤ 0, l ∈ I . N

Clearly,

ED

M

Note that the delayed argument µ(s|ϑ)−h ¯ depends on s and ϑ. Define a new function sdelay (s|ϑ), for each admissible duration vector ϑ ∈ V, n o sdelay (s|ϑ) = sup  ∈ [−h, 3] : µ(|ϑ) ¯ = µ(s|ϑ) ¯ − h , s ∈ [−h, 3]. (27) µ(s ¯ delay (s|ϑ)|ϑ) = µ(s|ϑ) ¯ − h.

(28)

xi,l (t − h) = xi,l (µ(s|ϑ) ¯ − h) = xi,l (µ(s ¯ delay (s|ϑ)|ϑ)) = x˜i,l (sdelay (s|ϑ)).

(29)

PT

Thus, for t ≥ 0, we obtain

AC

CE

Then, by substituting (3) and (28) into (26), we have  ˙i,l x˜ (s) = f˜i,l ( x˜i,l (s), x˜i,l (sdelay (s|ϑ)), ϑi , ζ i )       = ϑ(i f i,l ( x˜i,l (s), x˜i,l (sdelay (s|ϑ)), ζ i ), µ(s ¯ delay (s|ϑ)|ϑ) ∈ (τi−1 , τi ], s ∈ (i − 1, i], i ∈ I3 , l ∈ IN ,    l x , s = 0, i = 1, l ∈ IN ,  i,l  0  x˜ (s) =  i−1,l  x ˜ (s), s = i, i = 2, 3, l ∈ IN ,     x˜1,l (s) = φl (s), s ≤ 0, l ∈ I .

(30)

N

Let x˜l (·|ϑ, ζ) := ( x˜1,l (·|ϑ1 , ζ 1 ), x˜2,l (·|ϑ2 , ζ 2 ), x˜3,l (·|ϑ3 , ζ 3 ))T denote the solution of system (30) corresponding to the admissible pair (ϑ, ζ). Let σ(s|ϑ) denote the unique integer such that ϑσ(s|ϑ)+1 > 0 and µ(s|ϑ) ¯ −h∈

X h σ(s|ϑ) i=0

ϑi ,

σ(s|ϑ)+1 X i=0

 ϑi .

(31)

To solve (30), we need give the explicit formula for sdelay (s|ϑ). Before giving a explicit formula for sdelay (s|ϑ), we first note that sdelay (s|ϑ), s ∈ [−h, 3] and for any given s ∈ [0, 3] and ϑ, there exists a unique i ∈ I3 such that ϑi > 0 and t ∈ [µ(i ¯ − 1|ϑ), µ(i|ϑ)). ¯ Now, we give the explicit formula for sdelay (s|ϑ) as a theorem below: 9

ACCEPTED MANUSCRIPT

10

/ 00 (2019) 1–??

Theorem 3. Let ϑ ∈ V, for each s ∈ [−h, 3],

 µ(s|ϑ) ¯ − h, if µ(s|ϑ) ¯ − τ < 0,   X    bsc+1 i  ϑ (s − bsc) + h + ϑ sdelay (s|ϑ) =    i=σ(s|ϑ)+1    σ(s|ϑ) + , otherwise. ϑσ(s|ϑ)+1

(32)

sdelay (s|ϑ) = µ(s ¯ delay (s|ϑ)|ϑ) = µ(s|ϑ) ¯ − h.

CR IP T

Proof. Suppose first that µ(s|ϑ) ¯ − h < 0. Based on the definition of sdelay (s|ϑ), we have µ(s ¯ delay (s|ϑ)|ϑ) < 0. Therefore, combining this equation with (22) yields (33)

Suppose now that µ(s|ϑ) ¯ ∈ [0, t f ), from (31), we have σ(s) > 0, and sdelay (s|ϑ) ∈ [σ(s), σ(s) + 1). In other words, σ(s) = bsdelay (s|ϑ)c, and it follows from (22) and (27) that µ(s ¯ delay (s|ϑ)|ϑ) =

σ(s) X

`

σ(s|ϑ)+1

ϑ +ϑ

`=1

(sdelay (s|ϑ) − σ(s|ϑ)) =

bsc X `=1

ϑ` + ϑbsc+1 (s − bsc) − h.

(34)

sdelay (s|ϑ) = σ(s|ϑ) +

AN US

It follows from sdelay (s|ϑ) ≤ s that we have 0 ≤ σ(s|ϑ) ≤ bsc. Thus, (34) can be rearranged to yields X ϑbsc+1 (s − bsc) + h + ϑi i=σ(s|ϑ)+1

σ(s|ϑ)+1

ϑ

.

This completes the proof.  Thus, formulas (9a), (9b) and (10) become, respectively,

≤ 0, ∀s ∈ [i − 1, i], j ∈ I5 , i ∈ I3 , l ∈ IN .

(35a) (35b)

N X 5 X ( x X n ylj (µ(s ¯ k |ϑ)) o x˜lj (sk |ϑ, ζ) ) ˜lj (sk |ϑ, ζ) 2 α − 2E + P el x∗j x∗j x∗j 5N N el l=1 |Ω | l=1 j=1 s ∈Ω k

PT

e ζ) = H(ϑ,

:= x j∗ −

x˜i,lj (s|ϑi , ζ i )

ED

g˜ i,lj ( x˜i,l (s|ϑi , ζ i ))

M

g˜ i,lj ( x˜i,l (s|ϑi , ζ i )) := x˜i,lj (s|ϑi , ζ i ) − x∗j ≤ 0, ∀s ∈ [i − 1, i], j ∈ I5 , i ∈ I3 , l ∈ IN ,

( X N X N X X X 5 X 5 x n ylj11 (µ(s ˜lj11 (sk1 |ϑ, ζ) ¯ k1 |ϑ)) ylj22 (µ(s ¯ k2 |ϑ)) o x˜lj22 (sk2 |ϑ, ζ) 1−α 4 · Cov , P el 2 x∗j1 x∗j1 x∗j2 x∗j2 25N 2 ( N l1 =1 l2 =1 s ∈Ω e l s ∈Ω e l j1 =1 j2 =1 l=1 Ω ) k1

k2

CE

n ylj1 (µ(s x˜lj1 (sk1 |ϑ, ζ) ¯ k1 |ϑ)) −2 1 ∗ Cov 1 ∗ , x j1 x j1

(36)

AC

Problem EOP becomes

ylj22 (µ(s ¯ k2 |ϑ)) o n ylj22 (µ(s ¯ k2 |ϑ)) o) E . x∗j2 x∗j2

e ζ) over V × Z3 subject to the continuous state inequality constraints (35a) and ] Problem EOP: Minimize H(ϑ, (35b). Note that the new time-scaling transformation has replaced the variable switching times τ1 and τ2 in the original ] Problem EOP with conventional decision parameter vector ϑi , i ∈ I3 , in the equivalent Problem EOP. 4.2. Constraint transcription technique Continuous state inequality constraints (35a) and (35b) are often handled by constraint transcription technique ] will be briefly discussed as follows. (see [33] and [38]). The application of this technique to Problem EOP 10

ACCEPTED MANUSCRIPT

11

/ 00 (2019) 1–??

Constraints (35a) and (35b) are equivalently to Gl (ϑi , ζ i ) :=

10 Z X j=1

i

i−1

n o max 0, g˜ i,lj ( x˜i,l (s|ϑi , ζ i )) ds = 0, i ∈ I3 , l ∈ IN .

(37)

G˜ lε,γ (ϑi , ζ i ) :=

10 Z X j=1

i

i−1

CR IP T

However, Gl (ϑi , ζ i ), i ∈ I3 , l ∈ IN , are non-smooth in (ϑi , ζ i ). In addition, constraint (37) does not satisfy the constraint qualification. Thus it is not advisable to compute it as such numerically. Thus, constraint (37) is slacked off to the following inequality constraint in the computation ϕε (˜gi,lj ( x˜i,l (s|ϑi , ζ i )))ds − γ ≤ 0, i ∈ I3 , l ∈ IN ,

where ε > 0, γ > 0 and 0, h i,l i2 g˜ j ( x˜i,l (s|ϑi , ζ i )) + ε 4ε g˜ i,lj ( x˜i,l (s|ϑi , ζ i )),

if g˜ i,lj ( x˜i,l (s|ϑi , ζ i )) < −ε, ,

if − ε ≤ g˜ i,lj ( x˜i,l (s|ϑi , ζ i )) ≤ ε, if g˜ i,lj ( x˜i,l (s|ϑi , ζ i )) > ε.

AN US

         i,l i,l i i ϕε (˜g j ( x˜ (s|ϑ , ζ ))) =        

(38)

It is shown in [33] that for every ε > 0, ∃ a γ0 (ε) > 0, such that if (38) is satisfied for any γ, 0 < γ < γ0 (ε), then the original constraints (35a) and (35b) are also satisfied. ] Now, by virtue of the results obtained in [33, 38], we obtain a sequence of approximate problems of Problem EOP as given below.

M

e ζ) over V × Z3 subject to constraints (38). Problem EOP: Minimize H(ϑ,

4.3. Gradient formulae for constraint function

PT

ED

Suppose that we have reached a (ϑ, ζ) ∈ V × Z3 such that the inequality constraints (38) are not satisfied. Then, P P3 ˜ l i i a gradient-based optimization method will be used to minimize N l=1 i=1 G ε,γ (ϑ , ζ ) with respect to (ϑ, ζ). The new (ϑ, ζ) obtained will be such that the constraints (38) are satisfied. For this, we will need the gradient formulas of the constraint functions G˜ lε,γ (ϑi , ζ i ), l ∈ IN with respect to ϑi , i ∈ I3 and ζ i , i ∈ I3 . They are presented in the following theorems. Theorem 4. Given l ∈ IN , i ∈ I3 , ε > 0 and γ > 0, the gradient of the constraint function G˜ lε,γ (ϑi , ζ i ) with respect to ϑi , is given by

CE

∂G˜ lε,γ (ϑi , ζ i ) ∂ϑi

=

10 Z X j=1

i

i−1

∂ϕε (˜gi,lj ( x˜i,l (s|ϑi , ζ i ))) ∂˜gi,lj ( x˜i,l (s|ϑi , ζ i )) ∂˜gi,lj

∂ x˜i,l

ω ˜ i,l (s|ϑi , ζ i )ds,

(39)

AC

where ω ˜ i,l (s|ϑi , ζ i ) ∈ R5 , l ∈ IN , is the solution of the following auxiliary time-delay system:

 i,l i i ˜i,l i,l  ∂ f˜i,l ( x˜i,l (s), x˜i,l (sdelay (s|ϑ)), ϑi , ζ i )  i,l  ˙˜ i,l (s) = ∂ f ( x˜ (s), x˜ (si,ldelay (s|ϑ)), ϑ , ζ ) ω ˜ (s) + ω   ∂ϑi ∂ x˜ (s)    i,l i,l i,l i i h i,l i ˜  ∂ f ( x ˜ (s), x ˜ (s (s|ϑ)), ϑ , ζ ) ∂ x ˜ (s  delay delay (s|ϑ)) ∂sdelay  + ω ˜ i,l (sdelay (s|ϑ)) + , ∀s ∈ (i − 1, i], i i,l  ∂sdelay ∂ϑ  ∂x ˜ (sdelay )     1,l 5  ω ˜ (s) = 0 ∈ R , ∀s ≤ 0,     ω ˜ i−1,l (i − 1) ∈ R5 , i = 2, 3, l ∈ IN . ˜ i,l (i − 1) = ω 11

(40)

ACCEPTED MANUSCRIPT

12

/ 00 (2019) 1–??

Proof. Let ζ i be arbitrary but fixed, and let e be the unit vector in R. Then, ∂ x˜i,l (s) x˜i,l (s|(ϑi )ξ , ζ i ) − x˜i,l (s|ϑi , ζ i ) = lim , i ξ→0 ∂ϑ ξ

(41)

Step 1 : Preliminaries Define

CR IP T

where (ϑi )ξ = ϑi + ξe. i ξ Let µ(·) ¯ denote µ(·|ϑ ¯ i ), µ¯ ξ (·) denote µ(·|(ϑ ¯ ) ), sdelay (·) denote sdelay (·|ϑi ) and sξdelay (·) denote sdelay (·|(ϑi )ξ ) for simplicity. To proceed, we will prove this theorem as follows.

(F i,l )ξ (ς) := f˜i,l ( x˜i,l (ς), x˜i,l (sξdelay (ς)), (ϑi )ξ , ζ i ).

For each real number ξ ∈ R, let ( x˜i,l )ξ (·) denote the function x˜i,l (·|(ϑi )ξ , ζ i ) and x˜i,l (·) denote the function x˜i,l (·|ϑi , ζ i ). Then, from (30), for each ξ ∈ R, we have Z s Z s ( x˜i,l )ξ (s) = x˜i,l (i − 1) + (F i,l )ξ (ς)dς = x˜i,l (i − 1) + f˜i,l ( x˜i,l (ς), x˜i,l (sξdelay (ς)), (ϑi )ξ , ζ i )dς. (42) i−1

AN US

i−1

Define

(Γi,l )ξ (s) = ( x˜i,l )ξ (s) − x˜i,l (s) =

Z

s

i−1

h i (F i,l )ξ (ς) − F i,l (ς) dς.

M

It follows from the mean value theorem that we obtain Z 1 ( ˜i,l ξ i ∂( f˜ηi,l )ξ ) ∂( f˜ηi,l )ξ h i,l ξ ∂( fη ) i,l ξ i,l i,l ξ i,l (Γ ) (ς) + ( x˜d ) (ς) − x˜d (ς) + ξ dη, (F ) (ς) − F (ς) = ∂ϑi ∂ x˜i,l ∂ x˜di,l 0

(43)

(44)

ED

where ( f˜ηi,l )ξ := f˜i,l ( x˜i,l (ς) + ηΓξ (ς), η(( x˜di,l )ξ (ς) − x˜di,l (ς)), ϑi + ξe, ζ i ) and ( x˜di,l )ξ (ς) − x˜di,l (ς) = =

x˜di,l (sξdelay (ς)|(ϑi )ξ , ζ i ) − x˜di,l (sdelay (ς)|ϑi , ζ i )

x˜di,l (sξdelay (ς)|(ϑi )ξ , ζ i ) − x˜di,l (sdelay (ς)|(ϑi )ξ , ζ i ) + x˜di,l (sdelay (ς)|(ϑi )ξ , ζ i ) − x˜di,l (sdelay (ς)|ϑi , ζ i )

PT

=

x˜di,l (sξdelay (ς)|(ϑi )ξ , ζ i ) − x˜di,l (sdelay (ς)|(ϑi )ξ , ζ i ) + (Γi,l )ξ (sdelay (ς)).

(45)

AC

CE

n o The state set ( x˜i,l )ξ (·) : ξ ∈ [−1, 1] is equi-bounded on [0, 3] due to the fact that f˜i,l ( x˜i,l (ς), x˜i,l (sdelay (ς)), ϑi , ζ i ) satisfies the linear growth condition. Therefore, there exists a real number c1 > 0 such that for each ξ ∈ [−1, 1], ( x˜i,l )ξ (s) ∈ N5 (c1 ), s ∈ [0, 3],

where N5 (c1 ) denotes the closed ball in R5 of radius c1 centered at the origin. It should be noted that N5 (c1 ) is convex, therefore, for each ξ ∈ [−1, 1], x˜i,l (s) + η(Γi,l )ξ (s) ∈ N5 (c1 ), s ∈ [0, 3], η ∈ [0, 1].

In addition, it is effortless to understand that for each ξ ∈ [−1, 1], ϑi + ηξe ∈ N(c2 ), where c2 := |ϑi | + 1. 12

ACCEPTED MANUSCRIPT

13

/ 00 (2019) 1–??

∂( f˜ηi,l )ξ ∂( f˜ηi,l )ξ are continuous. Consequently, it is clear from the comand ∂ x˜i,l ∂ x˜di,l pactness of N5 (c1 ) and N(c2 ) together with the definitions of φ that there exists a real number c3 > 0 such that, for each ξ ∈ [−1, 1],

From (3) and (30), it follows that

i,l ∂( f˜η )ξ ≤ c3 , ∂ x˜i,l 5×5

i,l ∂( f˜η )ξ ≤ c3 , ∂ x˜di,l 5×5

i,l ∂( f˜η )ξ ≤ c3 , ∂ϑi 5

l,ξ ∂φη ≤ c3 , s ∈ [0, 3], η ∈ [0, 1], ∂µ¯ 5

(46)

Step 2 : The function (Γi,l )ξ (s) is of order ξ.

CR IP T

l i where φl.ξ ¯ + ηξe) − h), and | · | denotes the Euclidian norm. η := φ (µ(s|ϑ

Let ξ ∈ [−1, 1] be arbitrary. Let s p denote a time point in the new time horizon such that µ(s ¯ p |ϑ) = h.

CE

PT

ED

M

AN US

It is clear from inequalities (43) and (44) that (Γi,l )ξ (s) 5 Z s Z 1 ( ∂( f˜ηi,l )ξ i,l ξ ∂( f˜ηi,l )ξ h i,l ξ = (Γ ) (ς) + x˜d (sdelay (ς)|(ϑi )ξ , ζ i ) − x˜di,l (sdelay (ς)|(ϑi )ξ , ζ i ) i,l 0 0 ∂ x˜i,l ∂ x˜d i ∂( f˜ηi,l )ξ ) i,l ξ ξ dηdς , +(Γ ) (sdelay (ς)) + 5 ∂ϑi Z s Z 1 Z s Z 1 ∂( f˜i,l )ξ ∂( f˜ηi,l )ξ i,l ξ η i,l ξ ≤ (Γ ) (ς)dηdς + (Γ ) (sdelay (ς))dηdς 0 0 ∂ x˜i,l 5 0 0 ∂ x˜i,l 5 d Z s Z 1 Z s Z 1 h i,l ξ i ∂( f˜ηi,l )ξ i,l i ξ i i ξ i x˜d (sdelay (ς)|(ϑ ) , ζ ) − x˜d (sdelay (ς)|(ϑ ) , ζ ) dηdς + + ξdηdς 5 0 0 ∂ϑi 0 0 5 Z s Z 1 Z Z s p 1 ∂( f˜i,l )ξ p ∂( f˜ηi,l )ξ i,l ξ η i,l ξ + (Γ ) (ς)dηdς = (Γ ) (s (ς))dηdς delay 5 0 0 5 ∂ x˜i,l ∂ x˜di,l 0 0 Z s Z 1 Z s Z 1 p p h i,l ξ i ∂( f˜ηi,l )ξ i,l i ξ i i ξ i x˜d (sdelay (ς)|(ϑ ) , ζ ) − x˜d (sdelay (ς)|(ϑ ) , ζ ) dηdς + + ξdηdς i 0 5 ∂ϑ 0 0 0 5 Z s Z 1 Z s Z 1 ∂( f˜i,l )ξ ∂( f˜ηi,l )ξ i,l ξ η i,l ξ + (Γ ) (ς)dηdς (Γ ) (s (ς))dηdς + delay 5 s p 0 ∂ x˜i,l s p 0 ∂ x˜i,l 5 d Z s Z 1 Z s Z 1 h ξ i ∂( f˜ηi,l )ξ ξ ξ + x˜(sdelay (ς)|θ, δ ) − x˜(sdelay (ς)|θ, δ ) dηdς + ξdηdς , s ∈ [0, 3]. i sp 0 5 s p 0 ∂ϑ 5

AC

• Case 1: When sdelay < 0 namely, s ∈ [0, s p ], taking the norm of both sides of (43) and applying the definition of c3 yield Z s Z 1 ( ) ∂( f˜ηi,l )ξ i,l ξ ∂( f˜ηi,l )ξ ∂( f˜ηi,l )ξ l,ξ l (Γi,l )ξ (s) = (Γ ) (ς) + (φη (ς) − φη (ς)) + ξ dηdς , i i,l i,l 0 0 5 5 ∂ϑ ∂ x˜ ∂ x˜ d

where

l l i i ¯ + ξe) − h) − φl (µ(ς|ϑ ¯ ) − h) = ξ φl,ξ η (ς) − φη (ς) = φ (µ(ς|ϑ

13

i i ∂φl (µ(ς|ϑ ¯ + ηξe) − h) ∂µ(ς|ϑ ¯ + ηξe) , η ∈ [0, 1]. ∂µ¯ ∂ϑi

ACCEPTED MANUSCRIPT

14

/ 00 (2019) 1–??

Then, l ¯ i i + ηξe) − h) ∂µ(ς|ϑ ¯ + ηξe) ξ c3 , η ∈ [0, 1]. φξη (ς) − φη (ς) ≤ ξ ∂φ (µ(ς|ϑ ≤ 3 ∂µ¯ ∂ϑi

(47)

Hence, we give

Z s (Γi,l )ξ (s) ≤ c3 ξ + 3c2 ξ + (Γi,l )ξ (ς) dς, s ∈ [0, s p ]. c 3 3 5 5 0

CR IP T

It follows from Gronwall’s Lemma (see Theorem 2.8.6 in [33]) that (Γi,l )ξ (s) ≤ ξ (c3 + 3c2 ) exp(c3 ), s ∈ [0, s p ]. 3 5

(48)

AN US

• Case 2: When sdelay ≥ 0, namely, s ∈ [s p , 3], (Γi,l )ξ (s) 5 Z s Z 1 ( ∂( f˜ηi,l )ξ i,l ξ ∂( f˜ηi,l )ξ h i,l ξ = x˜d (sdelay (ς)|(ϑi )ξ , ζ i ) − x˜di,l (sdelay (ς)|(ϑi )ξ , ζ i ) (Γ ) (ς) + i,l 0 0 ∂ x˜i,l ∂ x˜d i ∂( f˜ηi,l )ξ ) i,l ξ +(Γ ) (sdelay (ς)) + ξ dηdς , i 5 ∂ϑ Z s Z 1 Z Z s 1 ∂( f˜i,l )ξ ∂( f˜ηi,l )ξ i,l ξ η i,l ξ 2 (Γ ) (ς)dηdς + ≤ ξ (c3 + 3c3 ) exp(c3 ) + (Γ ) (sdelay (ς))dηdς 5 s p 0 ∂ x˜i,l s p 0 ∂ x˜i,l 5 d Z s Z 1 Z s Z 1 ∂( f˜i,l )ξ i ∂( f˜ηi,l )ξ h i,l ξ η i,l i ξ i i ξ i + ξdηdς . x˜d (sdelay (ς)|(ϑ ) , ζ ) − x˜d (sdelay (ς)|(ϑ ) , ζ ) dηdς + s p 0 ∂ x˜i,l 5 5 s p 0 ∂ϑi d

ED

M

Applying the definition of sdelay and sdelay > 0 yields Z s Z 1 ξ Z s Z s ∂ f˜η i,l ξ i,l ξ (Γ ) (s (ς))dηdς ≤ c (Γ ) (s (ς)) dς ≤ c3 (Γi,l )ξ (ς) dς. 5 delay 3 delay 5 ˜d sp 0 ∂ x sp 0 5 It follows from the mean value theorem that

(49)

PT

i,l ξ,η ∂ x˜ (sdelay (ς)|ϑi + ηξe, ζ i ) ∂sξ,η delay ξ i,l i,l i ξ i i ξ i x˜ (s ξ ≤ 3c ≤ (ς)|(ϑ ) , ζ ) − x ˜ (s (ς)|(ϑ ) , ζ ) 5 3 ξ , d delay d delay ∂sdelay ∂ϑi

ξ where sξ,η delay (ς) = sdelay (ς) + η(sdelay (ς) − sdelay (ς)), η ∈ [0, 1].

AC

CE

By Gronwall’s Lemma (see Theorem 2.8.6 in [33]), Z s (Γi,l )ξ (s) ≤ ξ (c3 + 3c2 ) exp(c3 ) + c3 ξ + 3c2 ξ + 2c3 Γξ (ς) dς 3 3 5 5 0 2 ≤ ξ (c3 + 3c3 )(exp(c3 ) + 1) exp(2c3 ), s ∈ [s p , 3].

It follows immediately from (48) and (50) that (Γi,l )ξ (s) ≤ ξ (c3 + 3c2 )(exp(c3 ) + 1) exp(2c3 ), s ∈ [0, 3]. 3 5

(50)

(51)

This inequality holds whenever the magnitude of ξ is sufficiently small (less than or equal to one) because of the fact that ξ ∈ [−1, 1] was arbitrary. Therefore, the function (Γi,l )ξ (s) is of order ξ. In view of (45)-(47), (49) and (51), we have ( x˜i,l )ξ (ς) − x˜i,l (ς) ≤ x˜i,l (sξ (ς)|(ϑi )ξ , ζ i ) − x˜i,l (sdelay (ς)|(ϑi )ξ , ζ i ) + (Γi,l )ξ (sdelay (ς)) d d d delay d 5 5 5 2 ≤ 3c3 ξ + ξ (c3 + 3c3 )(exp(c3 ) + 1) exp(2c3 ). (52) 14

ACCEPTED MANUSCRIPT

15

/ 00 (2019) 1–??

Step 3 : The definition of ρ and its properties Let the function ρ : [−1, 0) ∪ (0, 1] be defined as follow: ) −1 Z 3 ( λ1,ξ (s) + λ2,ξ (s) + λ3,ξ (s) ds, ρ(ξ) = ξ 5 5 5

(53)

0

where

1

∂ x˜i,l

0

− λ2,ξ (s) =

∂ f˜i,l ( x˜i,l (s), x˜di,l (s), ϑi , ζ i ) ∂ x˜i,l

Z

1

0

− 3,ξ

λ (s)

=

( ˜i,l i,l ∂ f ( x˜ (s) + ηΓξ (s), η(( x˜di,l )ξ (s) − x˜di,l (s)), ϑi + ξe, ζ i ) )

CR IP T

Z

× (Γi,l )ξ (s)dη,

∂ f˜i,l ( x˜i,l (s), x˜di,l (s), ϑi , ζ i ) ∂ x˜di,l i,l

)

∂ x˜di,l

× (( x˜di,l )ξ (s) − x˜di,l (s))dη,

( ˜i,l ∂ f ( x˜ (s) + ηΓξ (s), η(( x˜di,l )ξ (s) − x˜di,l (s)), ϑi + ξe, ζ i ) ∂ϑi 0 ) ∂ f˜i,l ( x˜i,l (s), x˜di,l (s), ϑi , ζ i ) − ξdη. ∂ϑi

Z

1

≤ 2c3 (c3 +

3c23 )(exp(c3 )

!

+ 1) exp(2c3 ) ,

PT

(55)

(56)

(57)

! 2 ≤ 2c3 3c3 + (c3 + 3c3 )(exp(c3 ) + 1) exp(2c3 ) ,

(58)

≤ 2c3 .

(59)

ED

−1 ξ λ2,ξ (s) −1 ξ λ3,ξ (s)

M

It follows from (46) and (51)-(52) that −1 ξ λ1,ξ (s)

(54)

( ˜i,l i,l ∂ f ( x˜ (s) + ηΓξ (s), η(( x˜di,l )ξ (s) − x˜di,l (s)), ϑi + ξe, ζ i )

AN US

λ1,ξ (s) =

CE

From (51)-(52), it follows that

x˜di,l (s)

+

x˜i,l (s) + η(Γi,l )ξ (s) →

η(( x˜di,l )ξ (s)

−

x˜di,l (s))

→

x˜i,l (s), as ξ → 0, x˜di,l (s),

as ξ → 0,

(60) (61)

AC

uniformly with respect to s ∈ [0, 1] and η ∈ [0, 1]. Furthermore, it is evident that ϑi + ηξe → ϑi , as ξ → 0,

uniformly with respect to η ∈ [0, 1].

(62)

Since the convergence in (60) and (61) take place inside the ball N5 (c1 ) , the convergence in (62) takes place ∂( f˜ηi,l )ξ ∂( f˜ηi,l )ξ ∂( f˜ηi,l )ξ are uniformly continuous on the compact set N5 (c1 ) × , and inside of the ball N(c2 ), ∂ϑi ∂ x˜i,l ∂ x˜i,l d

15

ACCEPTED MANUSCRIPT

16

/ 00 (2019) 1–??

N5 (c1 ) × N(c2 ), ∂ x˜i,l i,l i,l ξ ˜ ∂ f ( x˜ (s) + ηΓ (s), η(( x˜di,l )ξ (s) − x˜di,l (s)), ϑi + ξe, ζ i ) ∂ x˜di,l

∂ f˜i,l ( x˜i,l (s) + ηΓξ (s), η(( x˜di,l )ξ (s) − x˜di,l (s)), ϑi + ξe, ζ i ) ∂ϑi i i ∂ x˜i,l (sξ,η delay (ς)|ϑ + ηξe, ζ ) ∂sdelay ∂sξ,η delay (s) ∂ϑi

∂ f˜i,l ( x˜i,l (s), x˜di,l (s), ϑi , ζ i )

→

∂ x˜i,l i,l i,l ˜ ∂ f ( x˜ (s), x˜di,l (s), ϑi , ζ i )

→

∂ x˜di,l

, as ξ → 0, , as ξ → 0,

→

∂ f˜i,l ( x˜i,l (s), x˜di,l (s), ϑi , ζ i ) , as ξ → 0, ∂ϑi

→

∂ x˜i,l (sdelay (ς)|ϑi , ζ i ) , as ξ → 0, ∂sdelay

→

∂sdelay (s) , as ξ → 0, ∂ϑi

CR IP T

∂ f˜i,l ( x˜i,l (s) + ηΓξ (s), η(( x˜di,l )ξ (s) − x˜di,l (s)), ϑi + ξe, ζ i )

AN US

ξ uniformly with respect to s ∈ [0, 3], η ∈ [0, 1], where the formula sξ,η delay (s) := sdelay (s) + η(sdelay (s) − sdelay (s)) is the corresponding delayed time of the control ϑi + ηξe. These results along with (51) imply that −1 −1 −1 ξ λ1,ξ (s) → 0, ξ λ2,ξ (s) → 0, ξ λ3,ξ (s) → 0, as ξ → 0, (63)

uniformly with respect to s ∈ [0, 3]. Based on (57)-(59) and (63), the conditions of Lebesgue Dominated convergence theorem (see Theorem 2.6.4 in [33]) hold. Therefore, ( ) Z 3 −1 1,ξ −1 2,ξ −1 3,ξ lim ρ(ξ) = lim ξ λ (s) + ξ λ (s) + ξ λ (s) ds = 0. (64) 5

0 ξ→0

5

5

M

ξ→0

Step 4 : The final step: Comparing ξ−1 (Γi,l )ξ (·) with ω ˜ i,l (·|ϑi , ζ i )

AC

CE

=

(Γi,l )ξ (s) = ( x˜i,l )ξ (s) − x˜i,l (s) Z s Z 1 ( ˜i,l ξ i ∂( f˜ηi,l )ξ ) ∂( fη ) i,l ξ ∂( f˜ηi,l )ξ h i,l ξ i,l (Γ ) (ς) + ( x ˜ ξ dηdς ) (ς) − x ˜ (ς) + d d ∂ϑi ∂ x˜i,l ∂ x˜di,l 0 0 ) Z s(X Z s ˜i,l i,l 3 ∂ f ( x˜ (ς), x˜di,l (s), ϑi , ζ i ) i,l ξ i,ξ (Γ ) (ς)dς λ (ς) dς + ∂ x˜i,l 0 0 i=1 Z s ˜i,l i,l Z s ˜i,l i,l i ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i ) h i,l ξ ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i ) i,l + ( x ˜ ξdς ) (ς) − x ˜ (ς) dς + d d ∂ϑi ∂ x˜di,l 0 0 ) Z s(X Z s ˜i,l i,l Z s ˜i,l i,l 3 ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i ) i,l ξ ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i ) i,ξ λ (ς) dς + (Γ ) (ς)dς + ξdς ∂ϑi ∂ x˜i,l 0 0 0 i=1 Z s ˜i,l i,l i ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i ) h i,l ξ + x˜d (sdelay (ς)|(ϑi )ξ , ζ i ) − x˜di,l (sdelay (ς)|(ϑi )ξ , ζ i ) + (Γi,l )ξ (sdelay (ς)) dς i,l ∂ x˜d 0 ) Z s(X Z s ˜i,l i,l Z s ˜i,l i,l 3 ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i ) ξ ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i ) i,ξ ξdς λ (ς) dς + Γ (ς)dς + ∂ϑi ∂ x˜i,l 0 0 0 i=1 ξ,η " # Z s ˜i,l i,l i i ∂ x˜i,l (sξ,η ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i ) i,l ξ delay (ς)|ϑ + ηξe, ζ ) ∂sdelay + (Γ ) (s (ς)) + ξ dς. (65) delay ∂sdelay ∂ϑi ∂ x˜di,l 0

PT

=

ED

Let ξ ∈ [−1, 0) ∪ (0, 1] be arbitrary but fixed. Next, it follows from (43)-(45), (54)-(56) and the mean value theorem that

=

=

16

ACCEPTED MANUSCRIPT

17

/ 00 (2019) 1–??

Integrating the auxiliary system (40) yields (where x˜di,l = x˜i,l (sdelay (ς|ϑ))) Z s ( ˜i,l i,l ∂ f ( x˜ (ς), x˜i,l (sdelay (ς|ϑ)), ϑi , ζ i ) i,l ∂ f˜i,l ( x˜i,l (ς), x˜i,l (sdelay (ς|ϑ)), ϑi , ζ i ) ω ˜ i,l (s|ϑi , ζ i ) = ω ˜ (ς) + ∂ϑi ∂ x˜i,l (s) 0 ) ∂ f˜i,l ( x˜i,l (ς), x˜i,l (sdelay (ς|ϑ)), ϑi , ζ i ) h i,l ∂ x˜i,l (sdelay (ς|ϑ)) ∂sdelay i + dς. (66) ω ˜ (sdelay (ς|ϑ)) + ∂sdelay ∂ϑi ∂ x˜i,l (sdelay )

0

+

+

Z

s

0

s

∂ f˜i,l ( x˜i,l (ς), x˜di,l (ς), ϑi , ζ i ) ∂ x˜i,l

0

i=1 ˜i,l i,l

" ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i )  ∂ x˜di,l

 ξ−1 (Γi,l )ξ (sdelay (ς)) − ω ˜ i,l (sdelay (ς)|ϑi , ζ i ))

ξ,η i i ∂ x˜i,l (sξ,η delay (ς)|ϑ + ηξe, ζ ) ∂sdelay

∂ϑi

∂sdelay

" # ξ−1 (Γi,l )ξ (ς) − ω ˜ i,l (ς|ϑi , ζ i ) dς

−

# ∂ x˜i,l (sdelay (ς)|ϑi , ζ i ) ∂sdelay (ς) dς. ∂sdelay ∂ϑi

AN US

=

ξ−1 (Γi,l )ξ (s) − ω ˜ i,l (s|ϑi , ζ i ) ) Z Z s (X 3 λi,ξ (ς) dς + ξ−1

CR IP T

Multiplying (65) by ξ−1 , and subtracting it from (66) yields

CE

PT

ED

M

Then, based on the definition of ρ(ξ) in (53), we have ξ−1 (Γi,l )ξ (s) − ω ˜ i,l (s|ϑi , ζ i ) 5 " # ) Z s ˜i,l i,l 3 Z s ( X ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i ) −1 i,l ξ i,l i i −1 i,ξ ξ (Γ ) (ς) − ω ˜ (ς|ϑ , ζ ) dς = ξ λ (ς) dς + ∂ x˜i,l 0 0 i=1 " Z s ˜i,l i,l  ∂ f ( x˜ (ς), x˜di,l (ς), ϑi , ζ i )  −1 i,l ξ + ξ (Γ ) (sdelay (ς)) − ω ˜ i,l (sdelay (ς)|ϑi , ζ i )) i,l ∂ x˜d 0 # i,l ξ,η ∂ x˜ (sdelay (ς)|ϑi + ηξe, ζ i ) ∂sξ,η ∂ x˜i,l (sdelay (ς)|ϑi , ζ i ) ∂sdelay (ς) delay − dς + 5 ∂sdelay ∂ϑi ∂sdelay ∂ϑi Z s Z s −1 i,l ξ −1 i,l ξ i,l i i i,l i i ˜ (ς|ϑ , ζ ) dς ≤ ρ(ξ) + ˜ (ς|ϑ , ζ ) dς + c3 ξ (Γ ) (ς) − ω c3 ξ (Γ ) (ς) − ω 5 5 0 0 Z s ∂ x˜i,l (sξ,η (ς)|ϑi + ηξe, ζ i ) ∂sξ,η ∂ x˜i,l (sdelay (ς)|ϑi , ζ i ) ∂sdelay (ς) delay delay + c3 − dς ∂sdelay ∂ϑi ∂sdelay ∂ϑi 5 0 Z s −1 i,l ξ i,l i i ≤ ρ(ξ) ˜ + 2c3 ξ (Γ ) (ς) − ω ˜ (ς|ϑ , ζ ) dς, 5 0 where

AC

ρ(ξ) ˜ := ρ(ξ) +

Clearly,

Z

0

s

i,l ξ,η ξ,η i i ∂ x˜ (sdelay (ς)|ϑ + ηξe, ζ ) ∂sdelay ∂ x˜i,l (sdelay (ς)|ϑi , ζ i ) ∂sdelay (ς) c3 − dς. ∂sdelay ∂ϑi ∂sdelay ∂ϑi 5 ρ(ξ) ˜ → 0, as ξ → 0.

It follows from Gronwall’s Lemma (see Theorem 2.8.6 in [33]) that ξ−1 (Γi,l )ξ (s) − ω ˜ i,l (s|ϑi , ζ i ) ≤ ρ(ξ) ˜ exp(2c3 ), s ∈ [0, 3]. 5

(67)

It should be noted that ξ ∈ [−1, 0) ∪ (0, 1] is arbitrary. Therefore, we can take the limit as ξ → 0 in (67) and apply (64) to get lim ξ−1 (Γi,l )ξ (s) = ω ˜ i,l (s|ϑi , ζ i ), s ∈ [0, 3].

ξ→0

17

(68)

ACCEPTED MANUSCRIPT

18

/ 00 (2019) 1–??

From (41), (43) and (68), it follows that ∂ x˜i,l (s) x˜i,l (s|(ϑi )ξ , ζ i ) − x˜i,l (s|ϑi , ζ i ) = lim = lim ξ−1 (Γi,l )ξ (s) = ω ˜ i,l (s|ϑi , ζ i ), s ∈ [0, 3]. i ξ→0 ξ→0 ∂ϑ ξ



(69)

Theorem 5. Given l ∈ IN , i ∈ I3 , ε > 0 and γ > 0, the gradient of the constraint function G˜ lε,γ (ϑi , ζ i ) with respect to ζ i , is given by

∂ζ i

:=

10 Z X

∂ϕε (˜gi,lj ( x˜i,l (s|ϑi , ζ i ))) ∂˜gi,lj ( x˜i,l (s|ϑi , ζ i ))

i

∂˜gi,lj

i−1

j=1

∂ x˜i,l

λ˜ i,l (s|ϑi , ζ i )ds,

CR IP T

∂G˜ lε,γ (ϑi , ζ i )

(70)

where λ˜ i,l (s|ϑi , ζ i ) ∈ R5×10 , i ∈ I3 , l ∈ IN , is the solution of the following auxiliary time-delay system:

(71)

AN US

 i,l i i i,l i i ˜i,l i,l ˜i,l i,l    ˜˙ i,l (s) = ∂ f ( x˜ (s), x˜ i,l (sdelay ), ϑ , ζ ) λ˜ i,l (s) + ∂ f ( x˜ (s), x˜ (sdelay ), ϑ , ζ ) λ˜ i,l (sdelay ) λ   i,l ∂ x˜ (s)   ∂x ˜ (sdelay )     ˜i,l ( x˜i,l (s), x˜i,l (sdelay ), ϑi , ζ i ) ∂ f  , ∀s ∈ (i − 1, i] +   ∂ζ i    1,l 5×10 ˜  λ (s) = 0 ∈ R , ∀s ≤ 0,     λ˜ i,l (i − 1) = λ˜ i−1,l (i − 1), i = 2, 3, l ∈ I . N

Proof: The proof is similar to that given for Theorem 3.  4.4. Fast approximation strategy

M

To calculate the gradients of the constraint function, it is clear from Theorems 3 and 4 that we need to calculate the solution of system (40) and the solution of system (71). This is computational expensive. Here, we propose a fast approximate scheme for the calculation the constraint functions and the gradients of the constraint functions. Let the time interval [i − 1, i], i ∈ I3 be equally divided into σi > 0, i ∈ I3 , subintervals, i.e., 1 = ~¯ i , i ∈ I3 , σi

ED

siki − siki −1 =

where ~¯ i is a constant; siki = i − 1 + ki ~¯ i , ki = 1, 2, . . . , σi , i ∈ I3 ; si0 = i − 1. Then, the constraint function G˜ lε,γ (ϑi , ζ i ) can be approximated as # ¯i 10 X σi " X ~ i,l i,l i i,l i,l i i i i i ϕε (˜g j ( x˜ (ski |ϑ , ζ ))) + ϕε (˜g j ( x˜ (ski −1 |ϑ , ζ ))) · − γ. 2 j=1 i

PT

G˜ lε,γ (ϑi , ζ i ) ≈

(72)

k =1

CE

For (72), it has been shown in [39] that for each ε > 0, if " # ¯i ~ ϕε (˜gi,lj ( x˜i,l (siki |ϑi , ζ i ))) + ϕε (˜gi,lj ( x˜i,l (siki −1 |ϑi , ζ i ))) · − ε/4 ≤ 0, ki ∈ Iσi , j ∈ I10 , i ∈ I3 , 2

(73)

AC

are satisfied, then the constraints (35a) and (35b) are satisfied at each discretization point. Similarly, the gradient formulas ∂G˜ lε,γ (ϑi , ζ i )/∂ϑi , i ∈ I3 , l ∈ IN , can be approximated as ∂G˜ lε,γ (ϑi , ζ i ) ∂ϑi +

≈

i,l i,l i 10 X σi " ∂ϕ (˜ X ˜ (ski |ϑi , ζ i ))) ∂˜gi,lj ( x˜i,l (siki |ϑi , ζ i )) ε g j (x j=1

ki =1

∂˜gi,lj

∂ x˜i,l

∂ϕε (˜gi,lj ( x˜i,l (siki −1 |ϑi , ζ i ))) ∂˜gi,lj ( x˜i,l (siki −1 |ϑi , ζ i )) ∂˜gi,lj

∂ x˜i,l

18

ξ˜i,l (siki |ϑi , ζ i )

# ¯i ˜ξi,l (si i |ϑi , ζ i ) · ~ . k −1 2

(74)

ACCEPTED MANUSCRIPT

19

/ 00 (2019) 1–??

Furthermore, the gradient formulas ∂G˜ lε,γ (ϑi , ζ i )/∂ζ i , i ∈ I3 , l ∈ IN , can be approximated as ∂G˜ lε,γ (ϑi , ζ i )

≈

∂ζ i +

i,l i,l i 10 X σi " ∂ϕ (˜ X ˜ (ski |ϑi , ζ i ))) ∂˜gi,lj ( x˜i,l (siki |ϑi , ζ i )) ε g j (x

∂ x˜i,l

∂˜gi,lj

j=1 ki =1

∂ϕε (˜gi,lj ( x˜i,l (siki −1 |ϑi , ζ i ))) ∂˜gi,lj ( x˜i,l (siki −1 |ϑi , ζ i )) ∂˜gi,lj

∂ x˜i,l

λ˜ i,l (siki |ϑi , ζ i )

# ¯i ˜λi,l (si i |ϑi , ζ i ) · ~ . k −1 2

(75)

Step 1 : Choose a (ϑ, ζ) = (ϑ1 , ϑ2 , ϑ3 , (ζ 1 )T , (ζ 2 )T , (ζ 3 )T )T from V × Z3 .

CR IP T

Based on (74)-(75), we propose an algorithm for computing the constraint functions and their gradients as follows. Algorithm 1

AN US

Step 2 : Solve system (30) using a fixed-step numerical integration scheme to obtain the solution x˜i,l (·|ϑi , ζ i ), i ∈ I3 , l ∈ IN . Step 3 : Solve systems (40) and (71) using a fixed-step numerical integration scheme to obtain the solution ξ˜i,l (·|ϑi , ζ i ) and λ˜ i,l (·|ϑi , ζ i ), i ∈ I3 , l ∈ IN , respectively. Step 4 : Calculate the constraint functions G˜ lε,γ (ϑi , ζ i ), i ∈ I3 , l ∈ IN , using (72). Step 5 : Compute the gradients of the constraint functions G˜ lε,γ (ϑi , ζ i ), i ∈ I3 , l ∈ IN , using (74) and (75).

Remark 3. The time horizon [i − 1, i], i ∈ I3 , are divided equally and a fixed-step numerical integration approach is used to solve the differential equations, where the partition nodes are directly regarded as integral nodes over a specific time horizon [i − 1, i], i ∈ I3 .

CE

PT

ED

M

4.5. Hybrid optimization algorithm Problem EOP is rather complex to allow for analytical solutions to be obtained. Thus, it is unavoidable to use numerical methods to solve this problem. In fact, Problem EOP is a smooth optimization problem and it could, in principle, be solved by using a gradient-based optimization technique (see, for example, [37]). However, these techniques are developed to search for local optimal solutions. Therefore, we combine gradient-based optimization methods with Genetic Algorithms (GAs) to develop a global search algorithm to solve Problem EOP. GAs [40] are adaptive heuristic search algorithm based on the evolutionary ideas of natural selection and genetics. As such, they represent an intelligent exploitation of a random search used to solve optimization problems. The traditional GAs are developed to handle unconstrained optimization problems. However, Problem EOP involves nonlinear continuous state inequality constraints on the state variables. We shall present a global optimization algorithm to solve Problem EOP, where a gradient-based optimization method is combined with a Genetic Algorithm. This algorithm is referred to as Algorithm 2. The parameters for Algorithm 2 are defined as follows: • σ ← (ϑ1 , ϑ2 , ϑ3 , (ζ 1 )T , (ζ 2 )T , (ζ 3 )T )T ∈ [v∗ , v∗ ] := V × Z3 . • nw is the number of genes.

AC

• i is the iteration index.

• w, w1 and w2 are the gene indexes.

• c¯ is the iteration index in crossover operator.

• M is the maximum number of iterations. • λ, rw1 and rw2 are random numbers taken from [0, 1]. • pc is the crossover probability. • nc is the total number of all offspring generated by crossover. 19

ACCEPTED MANUSCRIPT

20

/ 00 (2019) 1–??

• α˜ ∼ N(0, 1), β˜ ∼ N(0, 1). • pm is the mutation probability. • r is a random number. • nm is the total number of mutated offspring.

• % is the convergence tolerance of the cost function.

CR IP T

• M% is an integer for testing convergence (if the optimal objective value has not changed after M% iterations, then we terminate the algorithm).

Algorithm 2 Solve Problem EOP using gradient-based optimization and GAs.

Step 1 : Initialize data as follows: • Step 1.1 : Read known data such as %, M% , nm , pm , nc , nw , pc , M.

e l , xl , l ∈ N, h, z0 , W. • Step 1.2 : Read experimental data such as t f , D, N, Ωl , Ω 0

AN US

Step 2 : Set i ← 1 and randomly generate nw genes with a uniform distribution on V × Z3 , denoted by σiw .

Step 3 : For each w ∈ Inw , put σiw into Problem EOP, solve system (30) by using the Euler scheme together with Lagrange interpolation. Based on formula (72), check the values of G˜ lε,γ (σiw ), w ∈ Inw , l ∈ IN . If ∃w ∈ Inw and l ∈ IN such that G˜ lε,γ (σiw ) > 0 (i.e., σiw is infeasible for Problem EOP), then based on Algorithm 1 (an algorithm for computing the constraint functions and their gradients), minimizes G˜ lε,γ (σiw ) by using a gradientbased optimization algorithm to obtain a feasible point σ ˘ iw for Problem EOP (i.e., ∀l ∈ IN or w ∈ Inw , G˜ lε,γ (σ ˘ iw ) ≤ i i i 0) and set σw ← σ ˘ w . Calculate the value of the cost function (36), denoted as Jw , w ∈ Inw .

M

Step 4 : Crossover operator: • Step 4.1 : Set l˜ ← 0.

ED

• Step 4.2 : For each w1 ∈ Inw , if (rw1 < pc ) and ( for certain w2 ∈ Inw , w2 , w1 , rw2 < pc ), then

PT

l˜ ← l˜ + 1, σin

˜

w +l

l˜ ← l˜ + 1, σin

˜

w +l

← λσiw1 + (1 − λ)σiw2 ,

← (1 − λ)σiw1 + λσiw2 .

It is easy to check that the generated offsprings are still such that the inequality constraints (73) are satisfied.

AC

CE

• Step 4.3 : For each c¯ ∈ Inc , put σinw +¯c into Problem EOP, solve system (30) by using the Euler scheme together with Lagrange interpolation. Based on formula (72), check the values of G˜ lε,γ (σinw +¯c ), c¯ ∈ Inc , l ∈ IN . If ∃¯c ∈ Inc and l ∈ IN such that G˜ lε,γ (σinw +¯c ) > 0 (i.e., σinw +¯c is infeasible for Problem EOP), then based on Algorithm 1 (an algorithm for computing the constraint functions and their gradients), minimize G˜ lε,γ (σinw +¯c ) by using a gradientbased optimization algorithm to obtain a feasible point σ ˘ inw +¯c for Problem EOP (i.e., ∀l ∈ IN , G˜ lε,γ (σ ˘ inw +¯c ) ≤ 0) and set σinw +¯c ← σ ˘ inw +¯c . Calculate the value of the cost function (36), denoted as Jni w +¯c , c¯ ∈ Inc .

Step 5 : Mutation operator: • Step 5.1 : Inc is divided into two groups by σinw +∗ ← arg min Jni w +¯c , i.e. c¯ ∈Inc

Inc+ Inc−

n o ← c¯ ∈ Inc kσinw +l − σinw +∗ k ≥  , n o ← c¯ ∈ Inc kσinw +l − σinw +∗ k <  ,

where  is a small positive number and k · k represents the Euclidean norm. 20

ACCEPTED MANUSCRIPT

21

/ 00 (2019) 1–??

• Step 5.2 : Set m ˜ ← 0. • Step 5.3 : For each c¯ ∈ Inc , ( m ˜ ←m ˜ + 1, σinw +nc +m˜ ← ασ ˜ inw +¯c + (1 − α)σ ˜ inw +∗ , i i ˜ ˜ m ˜ ←m ˜ + 1, σnw +nc +m˜ ← βσnw +¯c + (1 − β)σinw +r ,

if rc¯ < pm and c¯ ∈ Inc+ , if rc¯ < pm and c¯ ∈ Inc− ,

where σis,nw +r is generated randomly from V × Z3 .

CR IP T

• Step 5.4 : For each m ∈ Inm , put σinw +nc +m into Problem EOP, solve system (30) by using the Euler scheme together with Lagrange interpolation. Based on formula (72), check the values of G˜ lε,γ (σinw +nc +m ), m ∈ Inm , l ∈ IN . If ∃m ∈ Inm and l ∈ IN such that G˜ lε,γ (σinw +nc +m ) > 0 (i.e., σinw +nc +m is infeasible for Problem EOP), then based on Algorithm 1 (an algorithm for computing the constraint functions and their gradients), minimize G˜ lε,γ (σinw +nc +m ) by using a gradient-based optimization algorithm to obtain a feasible point σ ˘ inw +nc +m for Problem l i i i EOP (i.e., ∀l ∈ IN , G˜ ε,γ (σ ˘ nw +nc +m ) ≤ 0) and set σnw +nc +m ← σ ˘ nw +nc +m . Calculate the value of the cost function i (36), denoted by Jnw +nc +m , m ∈ Inm .

AN US

Step 6 : Selection operator:

• Step 6.1 : Based on the value of Jdi , d ∈ Inw +nc +nm computed by Steps 3∼5, select the best nw genes from i nw + nc + nm ones, denoted by σi+1 w , w ∈ nw , where the smaller is the value of Jd , d ∈ Inw +nc +nm , the better are the particles. • Step 6.2 : Compute

Ji ←

min

d∈Inw +nc +nm

Jdi ,

σi ← arg

min

d∈Inw +nc +nm

Jdi .

M

    Step 7 : If i ≥ M or i ≥ M% and |J i − J i−M% +1 | < % , then output J i and σi , stop. Otherwise, set i ← i + 1 and go to Step 3.

ED

5. Numerical results

AC

CE

PT

Algorithms 1 and 2 are applied to seeking the optimal estimates of switching times and system parameters in Problem EOP and all computations were implemented in Matlab R2013a (The Mathworks Inc.) on Dell computer, which is equipped with Quad CPU (2 Core, 64-bit, clocked at 2.99 GHz) and 4 GB memory. Batch began at τ0 = 0 (hour) and ended at t f = 7.75 (hour). Euler algorithm together with Lagrange interpolation is used to solve system (3) with a step size of 1/3600 (hour). Some experimental data and known data are stated as follows: x01 := (0.102, 418.2609, 0, 0, 0)T , x02 := (0.2025, 441.337, 0, 0, 0)T , x03 := (0.173, 402.9348, 0, 0, 0)T , x04 := n o n o (0.2245, 509.8913, 0, 0, 0)T ; D := [2, 2.9] × [4, 5.8]; N := 4; Ω1 := 1, 2, 3, 4, 5, 6, 6.92 , Ω2 := 2, 3, 4, 5, 6, 6.5 , Ω3 := n o n o n o e 1 := 0.434, 0.867, 1.301, 1.734, 2.168, 2.601, 3 , Ω e 2 := 1, 2, 3, 4, 5, 6, 6.67 , Ω4 := 1.5, 2.5, 3.5, 4.75, 5.75, 6.75, 7.75 , Ω n o n o n e 3 := 0.450, 0.900, 1.349, 1.799, 2.249, 2.699, 3 , Ω e 4 := 0.581, 0.968, 1.355, 0.923, 1.385, 1.846, 2.308, 2.769, 3 , Ω o 1.839, 2.226, 2.613, 3 ; h := 0.26; Set the parameters values in Algorithm 2, M% := 50, % := 0.01, M := 600, pc := 0.8, pm := 0.2. To solve our robust optimization problem, we write a Matlab R2013a program based on Algorithms 1 and 2. Then, we solve our robust optimization problem through the application of our Matlab program for α = 0, 0.1, · · · , 1.0. For opt α = 0, 0.1, · · · , 1.0, the optimal estimates of the switching times τopt α and system parameters ζα are listed in Table opt opt 1. It is important to note that the system state with obtained optimal estimates (τα , ζα ) can fit to system state with different initial state vectors (x01 , x02 , x03 , x04 ).

21

ACCEPTED MANUSCRIPT

22

/ 00 (2019) 1–??

AN US

CR IP T

Table 1 The optimal estimates of the switching times and system parameters for α = 0, 0.1, · · · , 1.0. α ( τopt ζαopt α ) 0 ([2.974, 4.629]) [2.974,1.655,3.121,1.162,0.375,31.835,0.647,0.575,0.015,0.008,0.128,5.830,10.975,0.654,0.383,63.108,0.371, 0.709,0.019,0.008,0.132,6.679,16.669,0.527,0.354,38.039,0.635,0.996,0.015,0.010,0.058,8.128,14.457]T 0.1 ([2.516, 4.670]) [2.516,2.154,3.080,0.580,0.378,37.804,0.467,1.166,0.023,0.005,0.100,5.436,14.890,0.530,0.264,35.920, 0.326,1.403,0.009,0.010,0.067,8.407,11.762,0.682,0.299,45.178,0.532,0.785,0.012,0.005,0.104,7.881,13.883]T 0.2 ([2.551, 4.730]) [2.551,2.179,3.019,1.030,0.385,34.144,0.781,0.828,0.009,0.009,0.141,4.143,12.067,0.424,0.147,27.379, 0.385,0.614,0.021,0.007,0.147,4.267,5.904,1.046,0.314,50.072,0.493,1.128,0.017,0.006,0.088,6.451,7.191]T 0.3 ([2.729, 4.834]) [2.729,2.105,2.916,0.643,0.266,44.407,0.784,1.307,0.015,0.006,0.120,9.850,16.867,0.588,0.355,47.962, 0.557,0.853,0.021,0.007,0.117,4.110,12.580,1.106,0.141,44.020,0.422,1.057,0.009,0.007,0.134,8.655,15.969]T 0.4 ([2.578, 5.239]) [2.578,2.661,2.511,0.817,0.371,56.019,0.309,1.409,0.015,0.011,0.082,6.436,11.466,1.042,0.327,50.043, 0.599,1.163,0.014,0.007,0.145,5.913,12.111,0.587,0.327,27.347,0.620,1.186,0.018,0.010,0.150,5.627,14.320]T 0.5 ([2.790, 4.904]) [2.790,2.114,2.846,1.113,0.182,33.084,0.324,1.238,0.012,0.007,0.124,9.747,7.672,0.545,0.177,30.600, 0.454,1.431,0.013,0.005,0.117,6.350,9.796,1.170,0.400,35.458,0.597,0.869,0.016,0.011,0.095,4.549,17.383]T 0.6 ([2.660, 4.531]) [2.660,1.871,3.219,0.879,0.229,52.315,0.316,0.538,0.011,0.005,0.124,3.990,15.472,0.669,0.329,56.726, 0.607,0.956,0.013,0.011,0.079,5.368,7.222,0.918,0.166,25.119,0.411,0.824,0.016,0.010,0.108,5.927,13.557]T 0.7 ([2.891, 4.887]) [2.891,1.996,2.863,0.448,0.414,48.859,0.364,0.876,0.010,0.007,0.125,8.402,10.852,0.908,0.199,47.121, 0.556,1.126,0.011,0.011,0.126,6.900,10.478,0.501,0.212,42.256,0.787,0.781,0.022,0.005,0.053,9.547,15.528]T 0.8 ([2.665, 4.512]) [2.665,1.847,3.238,0.527,0.173,46.169,0.718,0.839,0.016,0.009,0.137,9.268,15.641,1.030,0.228,48.777, 0.304,0.493,0.019,0.005,0.144,9.330,12.861,0.626,0.239,32.699,0.502,0.625,0.019,0.012,0.128,7.549,7.627]T 0.9 ([2.474, 4.928]) [2.474,2.444,2.831,1.107,0.330,34.392,0.581,0.829,0.021,0.008,0.144,8.438,15.157,0.569,0.239,27.998, 0.490,1.151,0.015,0.008,0.085,6.124,14.086,0.600,0.160,67.651,0.679,1.096,0.010,0.005,0.086,9.715,9.219]T 1.0 ([2.981, 5.298]) [2.981,2.317,2.452,0.472,0.375,27.051,0.570,1.137,0.022,0.012,0.062,4.605,9.110,0.729,0.344,39.450, 0.430,1.292,0.008,0.004,0.051,4.681,9.575,1.027,0.143,66.013,0.641,0.531,0.013,0.008,0.118,3.933,11.356]T

AC

CE

PT

ED

M

For α = 0, 0.1, · · · , 1.0, and initial state vectors (x0l , l ∈ I4 ), to investigate solution robustness, we perturbed the optimal estimates for the switching times and system parameters σopt by a random percentage from the interval (0, 5%) in the positive and negative direction. We generate 100 random estimates according to the above procedure. e l , l ∈ I4 , generated For each perturbation, we generate 100,000 output realizations of the noisy measurement outputs Ω according to different distributions of output data. Then, we calculate the mean and variance of the cost function corresponding to each perturbation of the optimal estimates under different initial state vectors (x01 , x02 , x03 , x04 ). Our results are shown as box plots in Figs. 1-4. From Figs. 1-4, for different initial state vectors (x01 , x02 , x03 , x04 ), the tradeoff between error mean and error variance is clearly apparent. It is clear from Figs. 1-4 that with α increasing from 0 to 1, the error mean gradually drops while the error variance gradually increase. Imposing a small penalty on the error variance hardly changes the error mean, but significantly improves solution robustness because the deviation in the error variance is much more than the deviation in the error mean. This shows the solutions of our optimization problem are robust. It follows from Figs. 1-4 that we select α = 0.1 to compare with the previous results [7] to illustrative the effectiveness of the obtained results. This shows the advantages of our optimization problem: the solution robustness is distinctly on the increase at the expense of a insignificant cost to the relative error mean. For comparison, the relative errors between the computational values and the experimental data with different initial glycerol and biomass concentrations in this work and the ones in [7] are listed in Table 1, in which the relative errors are defined by X opt l |xlj (tk |τopt α , ζα ) − y j (tk )| elj :=

tk ∈Ωl

X

tk ∈Ωl

|ylj (tk )|

, l ∈ I4 , j ∈ I3 .

Table 1 shows the relative errors fall off significantly compared with the ones in [7]. Furthermore, the concentration opt changes computed by the optimal optimal estimates (τopt α , ζα ) is plotted in Fig. 5, in which the horizontal axes stand for time; the left vertical axes denote the concentrations of glycerol and 1,3-PD, while the right vertical axes stand for biomass; the scattergrams represent experimental results and the dashed lines display simulating curves. Fig. 5 shows 22

ACCEPTED MANUSCRIPT

23

/ 00 (2019) 1–??

Mean Variance

5

4

3

2

1

0 0 M

V

0.1 M

V

0.2 M

V

0.3 M

V

0.4 M

V

0.5 M

V

0.6 M

V

0.7 M

V

0.8 M

V

0.9 M

V

1.0 M

V

AN US

α

CR IP T

The values of mean or variance

6

Fig. 1. Mean (M, for short) and variance (V, for short) of the cost function for 100,000 output realizations generated according to different distribution with initial state vector x01 under 100 perturbations of the optimal estimates for the switching times and system parameters. (Please see the web version of this article to interpret the references to color.)

M

opt that the nonlinear switched time-delay with the optimal optimal estimates (τopt α , ζα ) can legitimately characterize the batch process with different initial conditions.

Table 1. The relative errors between this paper and literature [7].

ED

Substances j = 1 (Biomass) j = 2 (Glycerol) j = 3 (1,3-PD) j = 1 (Biomass) j = 2 (Glycerol) j = 3 (1,3-PD) j = 1 (Biomass) j = 2 (Glycerol) j = 3 (1,3-PD) j = 1 (Biomass) j = 2 (Glycerol) j = 3 (1,3-PD)

PT

l=1

l=3

l=4

AC

CE

l=2

elj (%) in the paper 7.53 4.35 5.18 11.44 6.55 12.29 18.39 6.48 12.42 12.97 5.94 17.44

elj (%) in [7] 41.87 10.21 19.11 38.14 13.31 27.87 43.25 12.02 15.94 23.47 16.26 24.16

6. Discussions and conclusions This paper studied a batch culture described by a nonlinear time-delay switched system, in which the output measurements are noisy due to system noise and measurement inaccuracies. The unknown switching times and unknown system parameters are to be determined. To identify the unknown quantities, we formulate the problem as a robust optimization problem governed by the switched time-delay system with continuous state inequality constraints and 23

ACCEPTED MANUSCRIPT

24

/ 00 (2019) 1–??

5

Mean Variance

The values of mean or variance

4.5 4

3.5

CR IP T

3

2.5 2

1.5 1 0.5 0 V

0.1 M

V

0.2 M

V

0.3 M

V

0.4 M

0.5 M

0.6 M

0.7 M

0.8 M

AN US

0 M

V

V

V

V

V

0.9 M

V

1.0 M

V

α

ED

M

Fig. 2. Mean (M, for short) and variance (V, for short) of the cost function for 100,000 output realizations generated according to different distribution with initial state vector x02 under 100 perturbations of the optimal estimates for the switching times and system parameters. (Please see the web version of this article to interpret the references to color.)

6

4

CE

3

PT

The values of mean or variance

5

Mean Variance

AC

2

1

0 0 M

V

0.1 M

V

0.2 M

V

0.3 M

V

0.4 M

V

0.5 M

V

0.6 M

V

0.7 M

V

0.8 M

V

0.9 M

V

1.0 M

V

α

Fig. 3. Mean (M, for short) and variance (V, for short) of the cost function for 100,000 output realizations generated according to different distribution with initial state vector x03 under 100 perturbations of the optimal estimates for the switching times and system parameters. (Please see the web version of this article to interpret the references to color.)

24

ACCEPTED MANUSCRIPT

25

/ 00 (2019) 1–??

4.5

Mean Variance

The values of mean or variance

4 3.5 3

2

1.5 1 0.5 0

0 M

V

0.1 M

V

0.2 M

V

0.3 M

V

0.4 M

V

0.5 M

V

0.6 M

V

0.7 M

V

0.8 M

V

0.9 M

V

1.0 M

V

AN US

α

CR IP T

2.5

Fig. 4. Mean (M, for short) and variance (V, for short) of the cost function for 100,000 output realizations generated according to different distribution with initial state vector x04 under 100 perturbations of the optimal estimates for the switching times and system parameters. (Please see the web version of this article to interpret the references to color.)

PT

7. Acknowledgements

ED

M

box constraints. Then, this problem is converted into a sequence of approximate subproblems. A global optimization algorithm is proposed to solve these approximate subproblems. Finally, numerical results show that optimal estimate of the switching times and systems parameters of the nonlinear switched time-delay system obtained when output measurements are carried in a noisy environment, can reasonably describe the actual process of batch culture. However, it should be note that the estimation method proposed in this paper is under assumption that the statistical properties of the output noise are available. Clearly, it is of practical significance to extend the estimation method to situation when the statistical properties of the measured output data are not available. This is an interesting further research project.

AC

CE

This work was supported by the Australian Research Council (Grant No. DP190103361), the National Natural Science Foundation of China (Grant Nos. 11771008, 11701063, 71831002 and 61773086), the Natural Science Foundation of Shandong Province in China (Grant Nos.: ZR2017MA005 and ZR2019MA031), the Fundamental Research Funds for the Centrol Universities (Grant Nos. DLMU3132019179, DLMU132019501 and DLMU3132019502 and DUT19LK37), the Project funded by China Postdoctoral Science Foundation (Grant No. 2019M651091), the Program for Innovative Research Team in University of Ministry of Education of China (Grant No. IRT-17R13), the Tianjin Thousand Talents Program and the Xinghai Project of Dalian Maritime University. References

[1] Biebl, H., Menzel, K., Zeng, A. P., Deckwer, W. D.: Microbial production of 1, 3-propanediol. Appl. Microbiol. Biot. 52(3), 289-297 (1999) [2] Zeng, A. P., Biebl, H.: Bulk chemicals from biotechnology: the case of 1, 3-propanediol production and the new trends. In Tools and Applications of Biochemical Engineering Science (pp. 239-259). Springer, Berlin Heidelberg, (2002) [3] Ashoori, A., Moshiri, B., Sedigh, A. K., Bakhtiari, M.R.: Optimal control of a nonlinear fed-batch fermentation process using model predictive approach. J. Process Contr. 19, 1162-1173 (2009) [4] Zheng, P., Wereath, K., Sun, J., van den Heuvel, J., Zeng, A. P.: Overexpression of genes of the dha regulon and its effects on cell growth, glycerol fermentation to 1, 3-propanediol and plasmid stability in Klebsiella pneumoniae. Process Biochem. 41(10), 2160-2169 (2006)

25

ACCEPTED MANUSCRIPT

26

2

100 0

0

2

4

M

4

ED

300 200 100 0

PT

glycerol and 1,3-PD(mmol/L)

c

400

0

0

6

t(h) 550 500

2

4

300 200

AN US

200

400

2

0

6

2

100 0

0

550 500

2

4

0

6

t(h) d

4

400 300

2

200 100 0

0

2

4

6

8

0

t(h)

CE

t(h)

4

AC

Fig. 5. (a) Comparisons of experimental data and simulating curves of biomass, glycerol and 1,3-PD concentrations with x01 ; (b) Comparisons of experimental data and simulating curves of biomass, glycerol and 1,3-PD concentrations with x02 ; (c) Comparisons of experimental data and simulating curves of biomass, glycerol and 1,3-PD concentrations with x03 ; (d) Comparisons of experimental data and simulating curves of biomass, glycerol and 1,3-PD concentrations with x04 .

26

biomass(g/L)

300

b

550 500

biomass(g/L)

400

CR IP T

4

biomass(g/L) glycerol and 1,3-PD(mmol/L)

a

550 500

biomass(g/L) glycerol and 1,3-PD(mmol/L)

glycerol and 1,3-PD(mmol/L)

/ 00 (2019) 1–??

ACCEPTED MANUSCRIPT

27

/ 00 (2019) 1–??

AC

CE

PT

ED

M

AN US

CR IP T

[5] Ye, J.X., Li, A., Zhai, J.G.: A measure of concentration robustness in a biochemical reaction network and its application on system identification. Appl. Math. Modell. 58, 270-280 (2018) [6] Rieckenberg, F., Ardao, I., Rujananon, R., Zeng, A. P.: Cell-free synthesis of 1, 3-propanediol from glycerol with a high yield. Eng. Life Sci. 14(4), 380-386 (2014) [7] Yuan, J.L., Zhu, X., Zhang, X., Yin, H.C., Feng, E.M., Xiu, Z.L.: Robust identification of enzymatic nonlinear dynamical systems for 1, 3-propanediol transport mechanisms in microbial batch culture. Appl. Math. Comput. 232, 150-163 (2014) [8] Cheng, G.M., Wang, L., Loxton, R.C., Lin, Q.: Robust optimal control of a microbial batch culture process. J. Optim. Theory. Appl. 16(1), 342-362 (2014) [9] Liu, C., Gong, Z, Lee, H.W.J., Teo, K.L.: Robust bi-objective optimal control of 1,3-propanediol microbial batch production process, J. Process Contr. 78, 170-182 (2019) [10] Wang, J., Ye, J.X. Yin, H.C. Feng, E.M., Wang, L.: Sensitivity analysis and identification of kinetic parameters in batch fermentation of glycerol. J. Comput. Appl. Math. 236(9), 2268-2276 (2012) [11] Wang, L., Feng, E.M., Xiu, Z.L.: Modeling nonlinear stochastic kinetic system and stochastic optimal control of microbial bioconversion process in batch culture. Nonlinear Anal. Model. 18(1), 99-111 (2013) [12] Zhang, J.X., Yuan, J.L., Feng, E.M., Yin, H.C., Xiu, Z.: Strong stability of a nonlinear multi-stage dynamic system in batch culture of glycerol bioconversion to 1,3-Propanediol, Math. Model. Anal., 21(2), 159-173 (2016) [13] An, Y., Tan, B., Wang, L., Chang, L., Wang, J.: Optimality condition and optimal control for a two-stage nonlinear dynamical system of microbial batch culture, Pac. J. Optim. 14(1), 1-13 (2018) [14] Zhai, J., Shen, B., Gao, J., Feng, E., Yin, H.: Optimal control of switched systems and its parallel optimization algorithm. J. Comput. Appl. Math. 261 , 287-298 (2014) [15] Yan, P., Ozbay, H.: Stability analysis of switched time delay systems. SIAM J. Control Optim. 47(2), 936-949 (2008) [16] Nestler, P., Sch¨oll, E., Tr¨oltzsch, F.: Optimization of nonlocal time-delayed feedback controllers. Comput. Optim. Appl. 64(1), 265-294 (2016) [17] G¨ollmann, L., Maurer, H.: Theory and applications of optimal control problems with multiple time-delays. J. Ind. Manag. Optim. 10, 413-441 (2014) [18] Lee, Y. S., Han, S.: An improved receding horizon control for time-delay systems. J. Optim. Theory Appl. 165(2), 627-638 (2015) [19] Liu, C., Loxton, R., Teo, K. L: A computational method for solving time-delay optimal control problems with free terminal time. Syst. Control Lett. 72, 53-60 (2014) [20] Denis-Vidal, L., Jauberthie, C., Joly-Blanchard, G.: Identifiability of a nonlinear delayed-differential aerospace model. IEEE T. Automat. Contr. 51, 154-158 (2006) [21] Ivanov, A. F., Mammadov, M. A., Trofimchuk, S. I.: Global stabilization in nonlinear discrete systems with time-delay. J. Glob. Optim. 56(2), 251-263 (2013) [22] Loxton, R., Teo, K.L., Rehbock, V.: An optimization approach to state-delay identification. IEEE T. Automat. Contr. 55, 2113-2119 (2010) [23] Silva, G.J., Datta, A., Bhattacharyya S.P.: PID controllers for time-delay systems. SIAM Rev. 47(4), 855 (2005) [24] Zamani, I., Shafiee, M.: Stability analysis of uncertain switched singular time-delay systems with discrete and distributed delays. Optimal Control Applications and Methods, 36(1), 1-28 (2015) [25] Liu, C.Y., Ryan, R., Lin, Q.,Teo, K.L.: Dynamic optimization for switched time-delay systems with state-dependent switching conditions, SIAM J. Control Optim. 56 (5), 3499-3523 (2018) [26] Liu, C.Y., Loxton, R., Teo, K.L.: Optimal parameter selection for nonlinear multistage systems with time-delays. Comput. Optim. Appl. 59(1-2): 285-306 (2014) [27] Liu, C.Y.: Modelling and parameter identification for a nonlinear time-delay system in microbial batch fermentation. Appl. Math. Model. 37, 6899-6908 (2013) [28] Yuan, J.L., Wang, L., Zhang, X., Feng, E.M., Yin, H.C., Xiu, Z.L.: Parameter identification for a nonlinear enzyme-catalytic dynamic system with time-delays, J. Glob. Optim. 62, 791-810 (2015) [29] Shao, J., Zhang, X., Feng, E.M., Xiu, Z.L.: Experimental design suboptimization for the enzyme-catalytic nonlinear time-delay system in microbial batch culture. J. Process Contr. 24, 1740-1746 (2014) [30] Yuan, J.L., Zhang, X., Liu, C.Y., Chang, L., Xie, J., Feng, E.M., Yin, H.C., Xiu, Z.L.: Robust optimization for nonlinear time-delay dynamical system of dha regulon with cost sensitivity constraint in batch culture, Commun. Nonlinear Sci. Numer. Simulat., 38, 140-171 (2016) [31] Lin, Q., Loxton, R., Xu, C., Teo, K. L. Parameter estimation for nonlinear time-delay systems with noisy output measurements. Automatica 60, 48-56 (2015) [32] Wu, D., Bai, Y., Yu, C.: A new computational approach for optimal control problems with multiple time-delay. Automatica 101, 388-395 (2019) [33] Teo, K.L., Goh, C.J., Wong, K.H.: A unified computational approach to optimal control problems. Long Scientific Technical, Essex, (1991) [34] Flury, B.: A First Course in Multivariate Statistics. New York: Springer, (1997) [35] Goberna, M.A., Lopez, ´ M,A.: Semi-infinite programming recent advances. Kluwer, Dordrecht, (2001) [36] Li, L., Yu, C., Zhang, N., Bai, Y.: A time-scaling technique for time-delay switched systems. Submitted to Discrete and Continuous Dynamical Systems-Series B [37] Lin, Q., Loxton, R., Teo, K. L.: The control parameterization method for nonlinear optimal control: a survey. J. Ind. Manag. Optim. 10(1), 275-309 (2014) [38] Loxton, R., Teo K.L., Rehbock, V., Yiu, K.F.C.: Optimal control problems with a continuous inequality constraint on the state and the control. Automatica 45, 2250-2257 (2009) [39] Martin, R., Teo, K.L.: Optimal control of drug administration in cancer chemotherapy. World Scientific, (1993) [40] Mitchell, M.: An introduction to genetic algorithms. MIT press, (1998) [41] Sun, Y. Q., Shen, J. T., Yan, L., Zhou, J. J., Jiang, L. L., Chen, Y., Feng, E.M., Yuan, J.L., Xiu, Z. L.: Advances in bioconversion of glycerol to 1, 3-propanediol: prospects and challenges. Process Biochem. 71, 134-146 (2018)

27