Input design for continuous time output error models. - Industrial

Subscriber access provided by NJIT | New Jersey Institute of Technology

Process Systems Engineering

Input design for continuous time output error models. Santhosh Kumar Varanasi, Chaitanya Manchikatla, and Phanindra Varma Jampana Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b05036 • Publication Date (Web): 01 Mar 2019 Downloaded from http://pubs.acs.org on March 2, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Input design for continuous time output error models Santhosh Kumar Varanasi, Chaitanya Manchikatla, and Phanindra Jampana∗ Department of Chemical Engineering, Indian Institute of Technology Hyderabad, Kandi, Sangareddy, Telangana- 502285, India. E-mail: [email protected] Phone: +91 9494424992 Abstract In this paper, the problem of input design for identification of continuous time output error models is considered. The input design problem is formulated as maximization of a measure of the Fisher Information Matrix, which defines the accuracy with which the system parameters can be estimated. The optimization problem involving the Fisher Information matrix depends on the true system parameters, which are unknown. Existing methods use an iterative approach to tackle this circular dependency. In this paper, the system transfer function is represented using generalized orthonormal basis functions which makes the Fisher Information matrix independent of the true system parameters. Therefore, the input design problem reduces to solving a single optimization problem instead of a series of optimizations that are performed in existing methods. Further, in existing approaches, the computational complexity of the optimization grows with the length of the input vector to be identified. In the current paper, the input is also represented using smooth basis functions, thereby making the complexity independent of the length of the input vector. A least squares approach for parameter estimation is given and it is shown that the estimates are unbiased for the

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

optimal inputs. The accuracy of the proposed method is demonstrated with the help of two examples. Further, input design for a continuous stirred tank reactor model and a quadruple tank experimental set up is considered in order to demonstrate the practical applicability of the new method.

1

Introduction

Mathematical models of physical systems are important in designing and tuning control algorithms. The models can be built from fundamental physical laws or from experimental data. Though physics based models give a complete description of the system, they are not easily obtained when there is a lack of a priori knowledge of the system. Furthermore, these models may have high computational demands thereby diminishing their appeal for real time control applications. System identification, which aims at identifying models from input-output data is often useful to overcome these difficulties 1 . In system identification, the covariance of the parameter estimates is an important measure of the accuracy of the algorithms. It is known that the covariance of the estimates obtained by Maximum Likelihood (ML) and Prediction Error (PE) methods is optimal. This optimal coherence is known as the Cramér-Rao lower bound (CRLB) and is equal to the inverse of the Fisher information matrix. The Fisher information matrix depends on the true (but unknown) parameters as well as the inputs provided to the system. To obtain parameter estimates with the smallest covariance, it is desirable to design the input in such a way that the inverse of the Fisher Information matrix becomes as small as possible. However, as the Fisher information matrix depends on the true parameters, the input design problem cannot be addressed directly 2 . In the existing literature, this issue has been dealt with in two different ways. The first approach is iterative input design where one iterates between the parameter estimation and input design 3,4 . The second approach is min-max or robust input design 5,6 carried out in the frequency domain. In the latter approach, a priori information, in terms of the parameters, which are assumed to lie in a 2


Page 2 of 31

Page 3 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


predefined region is required. In industrial processes, a limit on the amplitude of the input signal is often useful , e.g. mechanical limits on position of pneumatic valves and maximal velocity of hydraulic actuators result in such bounds. Input design in frequency domain with amplitude constraints is highly non-linear which makes the optimization problem computationally expensive 7 . The iterative input design is accomplished in the time domain, where amplitude constraints can be incorporated in a straight forward manner. Optimal input design for discrete time (DT) models using the iterative approach has been studied by various authors 7–9 . Even though discrete time models are ubiquitous in system identification literature, continuous time (CT) models 10–13 offer two main advantages : 1) CT models can handle nonuniformly sampled data and 2) CT models can effectively handle stiff systems i.e., systems with time constants of different orders of magnitude. In this paper, we present an approach for input design of CT systems by improving on the existing methods for DT systems 7 . The existing algorithms for DT systems suffer from three main issues. i) Dependency on true parameters of system, ii) Computational cost and iii) Need for multiple experiments. The computation of Fisher information matrix requires the parameters of transfer function which in general are unknown. To avoid this dependency, it is proposed to represent the unknown system transfer function using generalized orthonormal basis such as the TakenakaMalmquist functions 14 . By expanding the transfer function using this basis, the Fisher Information Matrix becomes independent of the true parameters of the system. Therefore, the input design reduces to solving a single optimization problem in contrast to the multiple optimizations that are performed in the iterative procedure. If the system contains dynamics which are separated by several orders of magnitude, the inputs also need to be sampled at a higher rate. Since the dimension of the optimization variable, i.e., the input vector, grows with faster sampling and the final time considered, existing methods can be computationally demanding. In order to handle this issue of higher computational complexity, the inputs can be represented using basis functions to make the

3



Page 4 of 31

input design problem independent of the sampling rate and the final time. In this paper, smooth basis functions are used to ensure that the method is also applicable to systems which contain input derivatives in their dynamics. The main contributions of the current paper are: 1) Identification of optimal input for CT output error models for multiple input multiple output (MIMO) systems by representing the transfer functions using generalized orthonormal basis and 2) Reducing computational cost by representing the input using smooth basis functions such as Fourier or B-splines. The proposed method relies only on the step test data of the plant for computing the optimal input. Therefore, the main advantage of the proposed method in practice is the very small experimental effort required for designing the optimal input. The rest of the paper is organized as follows. In Section 2, the problem formulation and an explanation of the generalized orthonormal basis are provided. The amplitude constrained optimal input design is explained in Section 3. Simulation results on two examples are shown in Section 4. Sections 5 and 6 present identification results on a simulated continuous stirred tank reactor (CSTR) model and quadruple tank experimental system respectively. The paper ends with concluding remarks in Section 7.

2 2.1

Preliminaries Problem Formulation

Consider a continuous linear time invariant system in output error form for MIMO case (with q inputs and m outputs),

xj (t) =

q X

Gi,j (p, βi,j )ui (t) j = 1, 2, . . . , m;

i=1

4


(1)

Page 5 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with measurements given by

yj (tk ) = xj (tk ) + ηj (tk ), k = 0, 1, 2, . . . , N

(2)

where Y N = {y1 (tk ), y2 (tk ), . . . , ym (tk ) k = 0, 1, 2, . . . , N } is the output data, p is the differential operator, Gi,j (p, βi,j ) i = 1, 2, . . . , q; j = 1, 2, . . . , m are unknown transfer functions, ηj (tk ) are i.i.d. zero mean Gaussian random variables and βi,j ∈ Rni.j is a vector which contains the parameters of (i, j)th transfer function Gi,j . The objective of optimal input design is to identify an amplitude constrained input |ui (t)|≤ Ci (t), i = 1, 2, . . . , q which maximizes the Fisher information matrix, M =E

T ! d d logf (Y N |β, u(τ ), τ ≤ tN ) logf (Y N |β, u(τ ), τ ≤ tN ) dβ dβ

(3)

where, β = {βi,j , i = 1, 2, · · · , q, j = 1, 2, · · · , m}, f (Y N |β, u(τ ), τ ≤ tN ) is the probability density function of Y N given the parameter vector β and the input sequence u(τ ), τ ≤ tN . The prediction error is j (tk ) = ηj (tk ) = yj (tk ) − yˆj (tk ) where yˆj (tk ) denotes the condiP tional expectation, E(yj (tk )|Y k−1 ) = qi=1 Gi,j (p, βi,j )ui (tk ). Note that the parameter vector β and the input history are required for computing the values yˆj (tk ), k = 0, 1, 2, · · · , N . We assume for simplicity that the errors (j (tk )) are i.i.d. zero mean Gaussian random variables with variance σ 2 . The Fisher information matrix can be then written as (proof attached in supporting information), " m N T # d 1 XX d yˆj (tk /β) yˆj (tk /β) Mβ = κ0 j=1 k=0 dβ dβ

(4)

where κ0 = σ 2 . Note that the method in this paper is applicable even when the errors do not have the same variance. However, in this case Mβ will be a weighted sum of the terms T d d y ˆ (t /β) y ˆ (t /β)) as shown in the supporting material. dβ j k dβ j k The computation of Fisher information matrix in Eq. (4) requires the partial derivatives 5



Page 6 of 31

of yˆj (tk ) at the true model parameters β. Since these are unknown, a sequential input design approach is generally used in which one iterates between estimating the parameters and the inputs. of the following steps. 1. A pseudo random binary sequence (PRBS) signal is considered as an input and the output data is obtained. 2. Using the input-output data, the parameters are first estimated and the optimal input is then obtained by maximizing the Fisher information for the estimated parameters. 3. Once the optimal input is identified, output data is generated with this input. This input-output data is again used in estimating the parameters. These steps are repeated until some termination criteria is met e.g. error between the estimated parameters in two consecutive iterations is in the tolerance range. In the iterative procedure, the input design is sensitive to the estimated system parameters. In order to avoid explicit dependence of parameters in the estimation of Fisher information matrix, it is proposed in this paper to expand the transfer function using generalized orthonormal basis such as the Takenaka-Malmquist functions.

2.2

Generalized orthonormal basis

The role of the basis function expansion is to make the input design problem independent of the parameters of the system transfer function. The k th order Takenaka-Malmquist basis functions in the Laplace domain are given by, √ k−1 2γk Y s − γi φk (s, γ¯ ) = s + γk i=0 s + γi

(5)

where γ¯ = (γ0 , γ1 , · · ·) is a sequence of poles satisfying γi ∈ R+ ∀ i = 1, 2, · · ·, to ensure stability of the identified system. Since Takenaka-Malmquist functions form a complete

6


Page 7 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


orthonormal basis of the Hardy space 14 , any stable transfer function can be expressed as

Gi,j (s, βi,j ) =

∞ X

l l αi,j φl (s, γ¯i,j )

(6)

l=0

l where, αi,j are the parameters of the system that need to be estimated. For practical pur-

poses, this sum is truncated at a finite number of terms. It may be noted that the truncation results in a set of functions that do not span the (infinite dimensional) Hardy space. These functions depend on the choice of the poles γi . If the dominant pole of the system is known and other poles are chosen so that they span the range of the desired frequencies, the truncation can still be effective in approximating the system. After truncation, the parameter L 1 0 ). It can be noted that the , · · · , αi,j , αi,j vector of the weights is represented as αi,j = (αi,j

weight vector αi,j depends on the true parameters βi,j .

3 3.1

Optimal input design and model identification Optimal input design

The predictions yˆ(tk ) can be obtained by applying Laplace transform on Eq. (2) with zero initial values. Representing the transfer function using truncated version of Eq. (6), the estimated model can be written as,

yˆj (tk ) =

q L X X

l αi,j ui,j,l f (tk )

(7)

i=1 l=0

l where ui,j,l ¯i,j ) ∗ ui (t) = f (t) = φl (t, γ

Rt 0

l l φl (t − τ, γ¯i,j )ui (τ )dτ where φl (tk , γ¯i,j ) represents the

l impulse response of filter φl (s, γ¯i,j ) and ∗ denotes the convolution operator. Therefore, the

Fisher information matrix is,

Mj (u) =

T N X ∂ yˆj (tk ) ∂ yˆj (tk ) k=0

∂α

7

∂α



Page 8 of 31

T

where α = αT1,1 αT2,1 . . . αTq,1 αT1,2 αT2,2 . . . αTq,2 . . . αT1,m αT2,m . . . αTq,m ∈ Rm∗q∗(L+1) , T l 0 2 L αi,j = αi,j is . The partial derivative of predictions with respect to αi,j αi,j . . . αi,j ∂ yˆj (tk ) = ui,j,l f (tk ) l ∂αi,j Z tk l φl (tk − τ, γ¯i,j )ui (τ )dτ = ≈

0 k X

l )ui (tr )δtr φl (tk − tr , γ¯i,j

r=0 l = uTi fi,j (k)

where ui (tr ), r = 0, 1, . . . , k is the sampled input data and for simplicity in notation we used l l l ui = (ui (t0 ), ui (t1 ), · · · , ui (tN ))T and fi,j (k) = (φl (tk − t0 , γ¯i,j ), φl (tk − t1 , γ¯i,j ), · · · , φl (tk − l tk , γ¯i,j ), 0, 0, · · · , 0)T ∈ RN +1 . Let Fi,j (k) be a N + 1 × L + 1 matrix whose lth column is given T ∂ yˆj (tk ) l by fi,j (k) for l = 0, 1, 2, · · · , L. Then, ∂αi,j = uTi Fi,j (k). Therefore,

∂ yˆj (tk ) ∂αj

T

= uT1 F1,j (k) uT2 F2,j (k) · · · , uTq Fq,j (k) T

where αj = α1,j , α2,j , · · · , αq,j

. Now, let i) Fj (k) be a q(N +1)×q(L+1) matrix such that T the blocks on the main diagonal are Fi,j (k) for i = 1, 2, · · · , q and ii) u = uT1 uT2 · · · , uTq . T ∂ yˆj1 (tk ) ∂ yˆj (tk ) It can then be seen that = uT Fj (k). Note that ∂α = 0 when j1 6= j2 . ∂αj j 2

Let Hj (k) be a q(N + 1) × qm(L + 1) block matrix such that the j th column block is Fj (k), i.e.,

Hj (k) = 0 0 · · · Fj (k) · · · 0 It can be noted that

∂ yˆj (tk ) ∂α

T

= uT Hj (k). Let Mj (u) =

8

∂ yˆj (tk ) k=0 ∂α

PN


∂ yˆj (tk ) ∂α

T

=

Page 9 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


PN

k=0

HTj (k)uuT Hj (k) for j = 1, 2, · · · , m. Then, the Fisher information matrix is given by

M(u) =

m X

Mj (u) =

j=1

m X N X

HTj (k)uuT Hj (k)

(8)

j=1 k=0

The objective is to find an input which maximizes some measure of conditioning of the Fisher information matrix. A widely used criteria known as D-optimality 7,8 is incorporated in this paper. 1

obj(u) = det[M(u)] q(N +1) := JD (M(u))

(9)

Remark. In the case when the initial conditions are non-zero, the predictor equation, Eq. (7) contains few more terms corresponding to the initial conditions. These terms depend on the parameters β, the initial states and initial inputs. Since these terms do not depend on the input dynamics, they will decay with time and do not contribute to the predictor equation for larger times. Hence, in order to account for the dependence of initial conditions in Eq. (1), the lower limit for k in Eq. (8) needs to be chosen to be a larger time value (instead of k = 0) depending on the effective time constant (τeff ) of system. Remark. If the overall duration of the process is relatively short, then the terms corresponding to the initial conditions may not be insignificant. Also, if obtaining the dependence on the initial conditions is an objective of the modelling process, then the method presented in this paper has limited applicability. This situation might arise in modelling of batch processes 15 . The method presented in this paper requires that the system is brought to a steady state before the optimal input is implemented on the system for parameter estimation. Remark. In this paper, we have only considered optimal input design for continuous time output error models. Maximizing the Fisher information matrix results in inputs that are optimal for the PE method as discussed earlier. In some applications, state space models, which are more descriptive in nature, are desired. For these models, parameter estimation using subspace methods is more common.

9



In discrete subspace identification, persistent excitation of inputs is sufficient for accurate identification (Lemma 10.4 and Lemma 10.5 in Verhagen et. al. 16 ). An advantage with the discrete subspace method is that the initial conditions do not influence input design. However, for continuous time state space models, conditions similar to persistent excitation on the input are not currently available. Therefore, optimal input design for continuous time subspace identification is more challenging.

3.2

Amplitude constrained identification problem

Amplitude constraints on the input signal can be incorporated into the optimization problem as, (P 1)

max JD (M(u)) u

subject to |ui (k)| ≤ Ci,k ∀i = 1, 2, . . . , q, k = 0, 1, · · · , N Let C be a diagonal matrix with elements Ci,k . Using the variable U = uuT , the optimization can be recast as, (P 2)

max JD (M(U)) u

subject to U ∈ Sn+ , U(i, j) ≤ D(i, j), rank(U) = 1 where Sn+ is the cone of positive semi-definite matrices and D is the matrix containing the appropriate products of amplitudes. Though the objective function in (P 2) is concave in the variable U, the feasible region is not convex due to presence of the rank constraint. In Manchester 7 , the rank constraint was removed and the authors obtained an input based on a randomized algorithm that guaranteed optimality of the Fisher Information Matrix (in the sense of Theorem. 1 given later) within a factor of π2 . In this method, the amplitude constraint is also dropped since it is enforced by the randomized technique.

10


Page 10 of 31

Page 11 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Instead of removing the rank constraint, a convex relaxation can also be considered. The rank constraint may be replaced by the nuclear norm (kUk∗ ), defined as the sum of singular values of a matrix. The nuclear norm has been known to promote low rank solutions 17 . It is equal to the trace of U when the matrix is positive semi-definite. However, it was noticed in the simulations that addition of this constraint did not result in a significant improvement. The objective function after removing the rank and the amplitude constraints is, (P 20 )

max JD (M(U)) U

subject to U ∈ Sn+ This problem is a maximization of a concave function over the convex cone of positive semi-definite matrices with linear constraints and can be solved efficiently. The solution to (P 20 ) is a (q ∗(N +1))×(q ∗(N +1)) matrix. In order to extract an input vector (ˆ u) of length q ∗ (N + 1), the following randomized procedure 7,18 is used. The first step is to compute a matrix D ≥ 0 such that U = DT D using singular value decomposition. The second step is to sample a vector ζ ∈ Rq∗(N +1) from the standard normal distribution. Then, the following solution vector is considered: ˆ = Csgn(DT ζ) u

(10)

where sgn is the sign or signum function. Theorem 1. The expectation of Fisher information matrix generated by Eq. (10) satisfies the following bound ˆ E(M(U))

2 M(U) π

(11)

ˆ that are at least û ˆ T = U) i.e., the randomized strategy is expected to give inputs (ˆ u | u 2/π times informative when compared with solution of problem (P 20 ) in terms of Fisher information matrix.

11



Page 12 of 31

Proof. An extension of the proof given in 7 is provided here. m X N X 2 2 T T ˆ − M(U) = ˆ ) − U Hj (k) E(M(U)) Hj (k) E(ˆ uu π π j=1 k=0

(12)

ˆ T ) − π2 U = CE(sgn(DT ζ)sgn(DT ζ)T )CT − π2 U = It can be noted that Q := E(ˆ uu CE(sgn(C−1 DT ζ)sgn(C−1 DT ζ)T )CT − π2 U. Using the fact that if x is Gaussian random variable and Σ is its covariance matrix, then E(sgn(x)sgn(x)T ) = observed that Q =

2 (Carcsin (C−1 UC−1 ) CT π

2 arcsin(Σ), π

it can be

− U). Now, since arcsin(X) X for any

X 0, we obtain that Q 0. Therefore, Q is positive semi-definite. ˆ Even though it was shown that E(M(U))

2 M(U) π

this does not necessarily imply

ˆ E(JD (M(U))) ≥ JD ( π2 M(U)). Therefore, it is not clear whether there exists an input of the ˆ = Csgn(DT ζ) which is optimal upto the form u

2 π

factor. However, it was observed in the

ˆ ≈ JD (M(U)) i.e., very ˆ for which JD (M(U)) numerical simulations that there exist inputs u close to being optimal. The main challenge with (P 20 ) is that the computational cost increases with the dimensionality of the matrix U. In order to reduce the complexity, the input (u) is exPNj T panded in an appropriate basis i.e., uj (tk ) = i=1 ai,j ψi,j (tk ) = aj ψj (tk ), where aj = T T a1,j a2,j . . . aNj ,j , ψj (tk ) = ψ1,j (tk ) ψ2,j (tk ) . . . ψNj ,j (tk ) . Hence, P T T uj = aj ψj (t0 ) ψj (t1 ) · · · , ψj (tN ) . Let Ψ be a q(N + 1) × q( qj=1 Nj ) block diagoT of size (N + 1) × Nj nal matrix with the diagonal blocks ψj (t0 ) ψj (t1 ) · · · , ψj (tN ) T T T T for j = 1, 2, · · · , q. It may be noted that u = Ψa where a = a1 a2 · · · aq and U = ΨAΨT with A = aaT . The parameter vector a contains the design variables and the objective of input design now reduces to estimating these variables.

12


Page 13 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The problem (P 20 ) can then be reposed as,

(P 3)

max JD (M(A)) A

subject to A ∈ Sn+ , where M(A) =

Pm PN j=1

k=0

ˆ is the solution to (P 3), then U ˆ = HTj (k)ΨAΨT Hj (k). If A

ˆ T . Decomposing U ˆ = DT D, we obtain the input u ˆ = Csgn(DT ζ) as before. ΨAΨ

3.3 3.3.1

Selection of poles and parameter estimation Pole selection

The input design proposed above requires knowledge of location of the poles for fixing the l . A method for selecting these poles is presented below. basis in Eq. (6), i.e., γ¯i,j

The poles are chosen based on the bandwidth of the system. To identify the bandwidth, sequential step tests are performed on the system i.e., in a multiple input case, we give step changes in one input and all the other inputs are held constant. This process is repeated for all the inputs. Since the step response plays a important role, instead of single step, a staircase signal as mentioned in 19 is used. The effective time constant is considered as the average of the estimated time constants from each of the step tests. From the effective time constant (τeff ) of the system, the bandwidth (ΩBW ) of system is considered as 1/time constant i.e., ΩBW = 1/τeff . The maximum frequency is considered as four times the estimated bandwidth and the poles are selected evenly in this frequency region 19 . Therefore, γmin = ΩBW and γmax = 4ΩBW and the poles are selected from a grid with these bounds and a pre-defined linear spacing.

3.3.2

Parameter estimation

Once the poles of the system are fixed and the input is identified, the parameters of the transfer function i.e., αj vector are obtained using a least squares approach. Using Eq. (7), 13



Page 14 of 31

a linear system of equations is obtained,

y = Aα + ϑ

(13)

where, T

y = y1 (t0 ) y1 (t1 ) . . . y1 (tN ) . . . ym (t0 ) ym (t1 ) . . .   A 0 . . . 0  1     0 A2 . . . 0    A= .  . . .  ..  . . . . . .     0 0 . . . Am  1,j,0 1,j,1 1,j,L q,j,0  uf (t0 ) uf (t0 ) . . . uf (t0 ) . . . uf (t0 )  1,j,0  u (t1 ) u1,j,1 (t1 ) . . . u1,j,L (t1 ) . . . uq,j,0 (t1 ) f f f  f Aj =  . . . .. . .  .. .. .. .. .. .   u1,j,0 (tN ) u1,j,1 (tN ) . . . u1,j,L (tN ) . . . uq,j,0 (tN ) f f f f ϑj = η1 (t0 ) η1 (t1 ) . . . η1 (tN ) . . . ηm (t0 ) ηm (t1 ) . . .

ym (tN )

uq,j,1 (t0 ) f

...



uq,j,L (t0 )  f

  uq,j,L (t ) 1  f  ..  .   q,j,1 q,j,L uf (tN ) . . . uf (tN ) T uq,j,1 (t1 ) . . . f .. .. . .

ηm (tN )

If AT A is invertible, the least squares solution is unique and is given by

ˆ = AT A α 

0   .. .   It may be noted that Mj (u) =  0   .. .  0

−1

0

˜ =u ˜u ˜ T is the solution to the optimization U

(14)

 ... 0  ..  .. . .   T . . . ATj Aj . . . 0  and M(u) = A A. If  .. ... . . . ..  . .  ... 0 ... 0 ˜ 0. problem (P 2) then it is clear that M(U)

0 ... .. . . . . 0 .. .

AT y

0 .. .

˜ is non-singular. Since the noise is assumed to be i.i.d. Therefore, the matrix AT A = M(U) 14


Page 15 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ˆ = α, which proves that the least squares zero mean Gaussian, it is easy to verify that E(α) estimator is unbiased. ˆ obtained using the randomized procedure, we know that E(M(U)) ˆ For the solution U 2 M(U) π

0, where the expectation is w.r.t. the random input realizations in Eq. (10).

Hence, E(AT A) is invertible. However, it is difficult to estimate the probability of invertibility of AT A. In the simulations performed in this paper, it was observed that the matrices ˆ were invertible. AT A obtained for all the realizations of u

4

Simulation Results

In this section, the proposed algorithm is implemented on two examples considered in Garnier et. al. 20 . The input is expanded using the Fourier and B-spline bases and a comparison between the two is given. In the Fourier Basis, since the major variations in frequency response occur in the range [γmin , γmax ], the input is chosen to contain many frequencies in this region. Two lower frequencies, i.e., ω = 0.001, 0.01 are also included in the input to obtain satisfactory identification in this regime as well. The details on construction of B-splines are given elsewhere 21,22 . In all the simulation results, the system is simulated in MATLAB with zero initial conditions to obtain noiseless {xj (tk ) for j = 1, 2, . . . , m and k = 0, 1, 2, . . . , N } data for a given input and Gaussian noise ηj (tk ) is added to xj (tk ) to obtain the output data {yj (tk ) for j = 1, 2, . . . , m and k = 0, 1, 2, . . . , N } so that the signal to noise ratio (SNR) is 20. SNR is defined as 10 × log10 (Psignal /Pnoise ) where Psignal , Pnoise are the power of the signal and noise respectively; power of a signal is defined as its root mean square value.

4.1

System-1

A two-input, one-output second order non-minimum phase system is considered.

15



x(t) =

−3p + 2 −0.5p + 1 u (t) + u2 (t) 1 p2 + 0.6p + 1 p2 + 4p + 3

A step test as mentioned in Section 3.3.1 is performed and the frequencies were computed to be γmin = [0.5429, 0.4343] for the two inputs. Six poles in the range [γmin , 4γmin ] for each transfer function are considered. After solving the problem (P 3) using the Fourier and BSpline basis for Ψ and obtaining the respective D matrices, different realizations of ζ are ˆ = Csgn(DT ζ) with Ci,k = 10 ∀i, k. A total of 500 used to compute the inputs using u realizations of the input sequences are obtained in this manner. For all the realizations, it was observed that the matrix AT A in Eq. (14) is invertible. The corresponding frequency responses of the estimated systems with Fourier and B-spline basis are shown in Figures 1 and 2 respectively. In these figures, the grey lines represent the identified systems with the different realizations of the inputs. 0

Magnitude (dB)

Magnitude (dB)

10 0 -10 -20

estimated system true system

270

180

90 10-2

-5 -10 -15 -20 -25 360

Phase (deg)

-30 360

Phase (deg)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 31

10-1

100

Frequency (rad/s)

180

90 10-2

101


270

10-1

100

Frequency (rad/s)

101

(a) Transfer function corresponding to Input-1 (b) Transfer function corresponding to Input-2

Figure 1: Comparison of frequency response for system-1 with Fourier basis. The grey lines represent the identified system with the different input and noise realizations. From Figures 1 and 2, it can be observed that there is a good match in frequency responses of the estimated and the true systems. However, with B-Splines basis it can be noted that the identification is not as accurate as with the Fourier basis at higher frequencies. In order to overcome this issue, the number of control points need to be increased. However, this 16


Page 17 of 31

0

Magnitude (dB)

Magnitude (dB)

10 0 -10 -20


270

180

90 10-2

-5 -10 -15 -20 -25 360

Phase (deg)

-30 360

Phase (deg)

10-1

100


270

180

90 10-2

101

Frequency (rad/s)

10-1

100

Frequency (rad/s)

101


Figure 2: Comparison of frequency response for system-1 with B-spline basis. The grey lines represent the identified system with the different input and noise realizations. results in a larger computational effort. Ratio of objective function to the upper bound

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 0.9 0.8 0.7 0.6 0.5 0.4 Fourier basis

B-spline basis

Figure 3: Box plot of ratio of objective function to the upper bound of Fourier and B-Spline basis for system-1 with different input realizations ˆ A box plot of ratio of objective functions for the various input realizations (i.e., M(U)) to the upper bound computed by using the optimal solution to (P 3) (i.e., M(U)) is shown in Fig 3. It can be observed from this figure that the median of ratio of objective function to the upper bound with Fourier and B-spline basis are 0.7206 and 0.8102 respectively. From this figure, it can also be noted that there exists an input in the Fourier basis whose objective value is very close to optimal (this is represented by the horizontal bar close to one in the box plot). The inputs for which the ratio is maximum are shown in Figures 4 and 5.

17



15

8

optimum input value

optimum input value

6 10

5

0

-5

4 2 0 -2 -4 -6 -8 -10

-10

-12 0

20

40

60

80

100

120

0

20

40

time

60

80

100

120

time

(a) Input-1

(b) Input-2

Figure 4: Identified optimal inputs for system-1 with Fourier basis. 10

optimum input value

10

optimum input value

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 31

5

0

-5

-10

5

0

-5

-10 0

20

40

60

80

100

120

0

20

time

40

60

80

100

120

time

(a) Input-1

(b) Input-2

Figure 5: Identified optimal inputs for system-1 with B-spline basis.

4.2

System-2

The second system considered is also a two-input one-output model,

x(t) =

−0.5p + 1 100 u (t) + u2 (t) 1 p2 + 0.6p + 1 p2 + 8p + 100

A step test as mentioned in Section 3.3.1 is performed and the frequenices were computed to be γmin = [0.5429, 3.6815]. Six and seven poles were used repectively for the two transfer functions in the range [γmin , 4γmin ]. The objective function value is computed by generating 500 input sequences using Eq. (10) and for all the realizations, it was observed that the matrix

18


Page 19 of 31

AT A in Eq. (14) is invertible. The corresponding frequency responses of the estimated systems with Fourier and B-spline basis are shown in Figures 6 and 7 respectively. 5

Magnitude (dB)

Magnitude (dB)

10

0

-10

Phase (deg)

Phase (deg)

0

-5

-10 360

-20 360 estimated system true system

270

180

90 10-2

-1

10 Frequency (rad/s)


180

0

-180 10-2

100

10-1

100

Frequency (rad/s)

101


Figure 6: Comparison of frequency response for system-2 with Fourier basis. The grey lines represent the identified system with the different input and noise realizations. 5

Magnitude (dB)

Magnitude (dB)

10

0

-10

Phase (deg)


270

180

90 10-2

0

-5

-10 360

-20 360

Phase (deg)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10-1

Frequency (rad/s)

0

-180 10-2

100


180

10-1

100

Frequency (rad/s)

101


Figure 7: Comparison of frequency response for system-2 with B-spline basis. The grey lines represent the identified system with the different input and noise realizations. From Figures 6 and 7, it can be noted that the estimated systems provide accurate frequency responses. A box plot of the ratio of objective functions for the various input realizations is shown in Fig 8. The median of the ratio of objective function to the upper bound with Fourier and B-spline basis are 0.7107 and 0.8335 respectively. Again, it may be noted that there exist input realizations that are close to being optimal. The inputs for which the ratio is maximum are shown in Figures 9 and 10. 19


Page 20 of 31

1 0.9 0.8 0.7 0.6 0.5 0.4 Fourier basis

B-spline basis

10

8

8

6

optimum input value

optimum input value

Figure 8: Box plot of ratio of objective function to the upper bound of Fourier and B-Spline basis for System-2 with different input realizations.

6 4 2 0 -2 -4 -6 -8

4 2 0 -2 -4 -6 -8 -10

-10

-12 0

20

40

60

80

100

120

0

20

40

time

60

80

100

120

100

120

time

(a) Input-1

(b) Input-2

Figure 9: Identified optimal inputs for system-2 with Fourier basis 10

10

8

8

optimum input value

optimum input value

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Ratio of objective function to the upper bound


6 4 2 0 -2 -4 -6

6 4 2 0 -2 -4 -6

-8

-8

-10

-10 0

20

40

60

80

100

120

0

20

time

40

60

80

time

(a) Input-1

(b) Input-2

Figure 10: Identified optimal inputs for system-2 with B-spline basis

5

CSTR Simulation Example

With a view towards application of the input design method for nonlinear processes around a steady state, a first principle based non-isothermal Continuous Stirred Tank Reactor (CSTR) 20 ACS Paragon Plus Environment

Page 21 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


model is considered.

5.1

Model description

The CSTR process consists of an irreversible, exothermic reaction, A → B, in a constant volume reactor cooled by a single coolant stream which is governed by the following equations:

q −E dCA = (CAf − CA ) − K0 CA exp dt V RT dT q(t) −∆H −E ρc Cpc −hA1 = (Tf − T ) + K0 CA exp + qc (t) 1 − exp (Tcf − T ) dt V ρCp RT ρCp V q c ρ c Cpc The dynamic model of CSTR was simulated in the MATLAB environment with the parameters 23 given in Table 1. The main objective of control is to maintain the concentration of reactant at a desired value with flow rate of coolant as the manipulated variable. The objective of input design is to identify a best possible input (coolant flow rate) from which the concentration predictions from the model can be made as accurate as possible. Table 1: Parameter values of CSTR Parameter

Value

Units

Process flow rate, q Feed concentration, CA0 Feed temperature, T0 Inlet coolant temperature, Tc0 CSTR volume, V Heat transfer term, hA Reaction rate constant, k0 Activation energy term, E/R Heat of reaction, δH Liquid densities, ρ, ρc Specific heats, Cp , Cpc

100 1 350 350 100 7 × 105 1 × 1011 1 × 104 −2 × 105 1 × 103 1

lit/min mol/lit K K lit cal/ min. K min−1 k cal/mol g/lit cal/g.K

21



5.2

Identification results

The proposed algorithm is applied on the CSTR model. To determine the bandwidth, firstly the system is brought to steady state and a step test (a two stair case signal each of magnitude 5 lit/min) is performed. The optimal inputs obtained from the proposed method with both Fourier and B-spline basis (γmin = 0.25) are shown in Figures 11(a) and 11(b) respectively. The approximate linear transfer function is then estimated using these identified optimal inputs. 8

15

6 10

optimum input value

4

optimum input value

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 31

2

0

-2

-4

5

0

-5

-6

-8

-10 0

5

10

15

20

25

30

0

time

5

10

15

20

25

30

time

(a) Fourier Basis

(b) B-Spline basis

Figure 11: Identified optimal inputs (coolant flow rate in lit/min) for CSTR system with Fourier basis and B-spline basis around a steady state coolant operating flow rate of 103 lit/min A maximum length (=7) pseudo random binary signal (PRBS) is then used to validate the identified models. The true and the estimated responses of the system are shown in Figure 12. From Figure 12, it can be observed that the identified models with both the B-splines and Fourier basis are able to predict the output response accurately (with root mean square errors (RMSE) of 0.0364 and 0.0442 respectively). The model obtained from the iterative procedure with a total of fifty iterations results in relatively large errors in the response with an RMSE of 0.1882. In all the methods, a third order approximation of the nonlinear model around the steady state have been considered.

22


Page 23 of 31

0.11 0.108 0.106

Concentration (mol/lit)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.104 0.102 0.1 0.098 0.096 0.094

response of Non-linear system Fourier basis model response B-spline basis model response iterative approach

0.092 0.09 0

10

20

30

40

50

time

Figure 12: Comparison of output response of true and estimated models with Fourier basis, B-spline basis and iterative approach for the CSTR model

6

Quadruple Tank experimental case study

In order to demonstrate the practical significance of the proposed method, input design on a quadruple tank experimental setup as shown in Figures 13(a) and 13(b) is considered.

6.1

Process description

In this experimental set up, the control valve openings (v1 and v2 ) are considered as inputs and the outputs are y1 and y2 , the water levels of Tanks 1 and 2 respectively. The governing equations of the system from mass balances are as follows 24,25

dh1 dt dh2 dt dh3 dt dh4 dt

a1 p 2gh1 + A1 a2 p =− 2gh2 + A2 a3 p =− 2gh3 + A3 a4 p =− 2gh4 + A4 =−

a3 p 2gh3 + A1 a4 p 2gh4 + A2 1 − γ2 f2 A3 1 − γ1 f1 A4

23


γ1 f1 A1 γ2 f2 A2

(15)


(a) Experimental setup of quadruple tank system

Page 24 of 31

(b) quadruple tank system flow sheet

Figure 13: Experimental Set-up and Flow sheet of quadruple tanks experiment where hi , ai and Ai are the level of liquid, outlet cross sectional area and cross sectional area of Tank i respectively. The outlet flow rate from the control valves 1 and 2 are denoted as f1 and f2 respectively. The parameters [γ1 , γ2 ] ∈ [0, 1] determine the fraction of flow that is received by the individual tanks. The measured output values are kc h1 and kc h2 . Linear control valves were used and hence the flow rate fi and the control valve opening vi are related by fi (t) = ki vi (t). The acceleration due to gravity is denoted by g and the parameters of the experimental set up are given in Table 2. The experimental operating points (as shown in Table 3) are designed in such a way that the system exhibits nonminimum phase characteristics. Table 2: Parameters of quadruple tank experimental set up Parameter

Symbol

unit

value

Cross sectional area of tank i Outlet cross sectional area of tank i Level sensor calibration constant Gravitational constant

Ai ai kc g

cm2 cm2 cm/sec2

176 2.55 2 981

24


Page 25 of 31

Table 3: Operating points of quadruple tank experimental set up Parameter

Symbol

Steady state heights h01 , h02 , h03 , h04 Steady state flow rates f r10 , f r20 Ratio of flow in Tank 1 to flow in Tank 4 γ1 Ratio of flow in Tank 2 to flow in Tank 3 γ2

6.2

unit

Operating point value

cm cm3 /sec -

10.1,7.9,4.6,2.5 284,350 0.43 0.34

Identification results

In order to determine the bandwidth of the system,  a series of step  tests were performed 0.0769 0.0357 and the frequencies are determined to be γmin =  . The number of poles 0.05 0.0357 considered for the four transfer functions were 3, 4, 4, 3 respectively. The identified optimal inputs with both Fourier and B-spline basis are shown in Figures 14 and 15 respectively. The transfer function is then estimated using these identified optimal inputs. 6

6

optimum input value

4

optimum input value

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2

0

-2

-4

4

2

0

-2

-4 -6

0

500

1000

1500

2000

2500

3000

0

500

time

1000

1500

2000

2500

3000

time

(a) Input-1

(b) Input-2

Figure 14: Identified optimal inputs for tanks in series experimental set up with Fourier basis The experiment is repeated with a maximum length (=11) PRBS signal and the transfer function is estimated using the iterative procedure as explained in Section 2. A multifrequency sinusoidal signal (with frequencies based on identified bandwidth) is then used to validate the identified models. The true and the estimated responses of the system are shown in Figure 16. 25


6

4

4

optimum input value

6

2

0

-2

-4

-6

Page 26 of 31

2

0

-2

-4

-6 0

500

1000

1500

2000

2500

3000

0

500

1000

time

1500

2000

2500

3000

time

(a) Input-1

(b) Input-2

Figure 15: Identified optimal inputs for tanks in series experimental set up with B-spline basis true output Fourier model B-spline model iterative approach

1

true output Fourier model B-spline model iterative approach

1

level in tank-2 in cm

1.5

level in tank-1 in cm

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

optimum input value


0.5 0 -0.5 -1

0.5

0

-0.5

-1 -1.5

0

500

1000

1500

2000

2500

3000

0

500

time in seconds

1000

1500

2000

2500

3000

time in seconds

(a) Response of Output-1

(b) Response of Output-2

Figure 16: Comparison of output response of true and estimated system with Fourier basis, B-spline basis and iterative approach for quadruple tank system From Figure 16, it can be observed that the identified model with B-splines is more accurate when compared to the model obtained using the Fourier basis. However there still exist minor errors (approximately 0.5 cm) at few time instants in response-1 for both the basis. The model obtained from the iterative procedure results in relatively larger errors in the responses. Input design using the iterative procedure is dependent on the orders of the system considered, which in general are unknown. In the current experimental study, the orders of the four transfer functions are fixed as (0,1),(0,2),(0,2) and (0,1) where the two numbers 26


Page 27 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


represent the input and output orders respectively. These orders are selected based on the linearized dynamics of the quadruple tank system. The large errors in the iterative method may be attributed to i) lower order models and ii) insufficient number of iterations performed for convergence. Since the iterative procedure involves repeated testing of the plant using the input designed in the previous iteration, the process is time consuming. Experimentally, a total of thirty six hours were spent performing the six iterations that were reported here. However, for the proposed method only one step test data was sufficient for designing the optimal input.

7

Conclusions

In the current paper, a new algorithm for amplitude constrained input design for continuous time output error models is presented. The optimal input design problem is formulated as maximization of the Fisher Information Matrix. This problem is non-convex and the optimization is dependent on the true parameters of system, which are unknown. The proposed method is based on a convex relaxation where the system transfer function is represented using a generalized orthonormal basis. This formulation makes the problem independent of true parameters of system. The only requirement is the selection of poles of the basis functions and a method of selecting the poles is provided. To reduce the computational complexity of the method, the input is also expanded using basis functions. A randomized strategy for obtaining the solution is presented and it is shown theoretically that these inputs are at least 2/π times as informative as the optimal inputs in the sense of expectation. A least squares approach for parameter estimation is given and it is shown that estimates are unbiased for the optimal inputs. Two examples were presented that demonstrated the accuracy of the frequency response obtained from the proposed method. It was also observed that the randomized strategy was able to provide inputs which are very close to being optimal. The method described in this paper is also applied for input design for a simulated CSTR

27



model and a quadruple tank experimental setup. The results show that the new method results in more accurate models compared to those obtained using the existing iterative input design method.

Supporting Information The derivation of Fisher information matrix for the multiple input multiple output (MIMO) case is presented.

References (1) Ljung, L. System identification: Theory for the user ; Prentice Hall, New Jersey, 1999. (2) Ng, T. S.; Goodwin, G. C.; Söderström, T. Optimal experiment design for linear systems with input-output constraints. Automatica 1977, 13, 571–577. (3) Barenthin, M.; Jansson, H.; Hjalmarsson, H. APPLICATIONS OF MIXED H2 AND H∞ INPUT DESIGN IN IDENTIFICATION. IFAC Proceedings Volumes 2005, 38, 458–463. (4) Gerencsér, L.; Hjalmarsson, H.; M˚ artensson, J. Identification of ARX systems with non-stationary inputsasymptotic analysis with application to adaptive input design. Automatica 2009, 45, 623–633. (5) Rojas, C. R.; Welsh, J. S.; Goodwin, G. C.; Feuer, A. Robust optimal experiment design for system identification. Automatica 2007, 43, 993–1008. (6) Rojas, C. R.; Aguero, J.-C.; Welsh, J. S.; Goodwin, G. C.; Feuer, A. Robustness in experiment design. IEEE Transactions on Automatic Control 2012, 57, 860–874. (7) Manchester, I. R. Input design for system identification via convex relaxation. 49th IEEE Conference on Decision and Control (CDC). 2010; pp 2041–2046. 28


Page 28 of 31

Page 29 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(8) Stojanovic, V.; Nedic, N.; Prsic, D.; Dubonjic, L. Optimal experiment design for identification of ARX models with constrained output in non-Gaussian noise. Applied Mathematical Modelling 2016, 40, 6676–6689. (9) Stojanovic, V.; Nedic, N. Robust identification of OE model with constrained output using optimal input design. Journal of the Franklin Institute 2016, 353, 576–593. (10) Garnier, H.; Wang, L. Identification of continuous-time models from sampled data; Springer, 2008. (11) Young, P. Parameter estimation for continuous-time modelsa survey. Automatica 1981, 17, 23–39. (12) Unbehauen, H.; Rao, G. Continuous-time approaches to system identificationa survey. Automatica 1990, 26, 23–35. (13) Rao, G. P.; Unbehauen, H. Identification of continuous-time systems. IEE ProceedingsControl theory and applications 2006, 153, 185–220. (14) Mi, W.; Qian, T.; Wan, F. A fast adaptive model reduction method based on Takenaka– Malmquist systems. Systems & Control Letters 2012, 61, 223–230. (15) Corbett, B.; Mhaskar, P. Subspace identification for data-driven modeling and quality control of batch processes. AIChE Journal 2016, 62, 1581–1601. (16) Verhaegen, M.; Verdult, V. Filtering and system identification: a least squares approach; Cambridge university press, 2007. (17) Candès, E. J.; Recht, B. Exact matrix completion via convex optimization. Foundations of Computational mathematics 2009, 9, 717. (18) Luo, Z.-Q.; Ma, W.-K.; So, A. M.-C.; Ye, Y.; Zhang, S. Semidefinite relaxation of quadratic optimization problems. IEEE Signal Processing Magazine 2010, 27, 20–34. 29



(19) Tangirala, A. K. Principles of system identification: Theory and practice; Crc Press, 2014. (20) Garnier, H.; Gilson, M.; Young, P. C.; Huselstein, E. An optimal IV technique for identifying continuous-time transfer function model of multiple input systems. Control engineering practice 2007, 15, 471–486. (21) Poyton, A.; Varziri, M. S.; McAuley, K. B.; McLellan, P.; Ramsay, J. O. Parameter estimation in continuous-time dynamic models using principal differential analysis. Computers & chemical engineering 2006, 30, 698–708. (22) Varanasi, S. K.; Jampana, P. Parameter Estimation and Model Order Identification of LTI Systems. IFAC-PapersOnLine 2016, 49, 1002–1007. (23) Morningred, J. D.; Paden, B. E.; Seborg, D. E.; Mellichamp, D. A. An adaptive nonlinear predictive controller. Chemical Engineering Science 1992, 47, 755–762. (24) Johansson, K. H. The quadruple-tank process: A multivariable laboratory process with an adjustable zero. IEEE Transactions on control systems technology 2000, 8, 456–465. (25) Vadigepalli, R.; Gatzke, E. P.; Doyle, F. J. Robust control of a multivariable experimental four-tank system. Industrial & engineering chemistry research 2001, 40, 1916–1927.

30


Page 30 of 31

Page 31 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical Abstract:

31


Input design for continuous time output error models. - Industrial

Recommend Documents