Economically Optimal Input Design Approach for Identification of

Subscriber access provided by Kaohsiung Medical University

Process Systems Engineering

Economically optimal input design approach for identification of constrained processes Abhishankar Kumar, Nabil Magbool Jan, and Sridharakumar Narasimhan Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b05187 • Publication Date (Web): 02 May 2018 Downloaded from http://pubs.acs.org on May 4, 2018

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 41 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

Economically optimal input design approach for identification of constrained processes Abhishankar Kumar,† Nabil Magbool Jan,†,‡ and Sridharakumar Narasimhan∗,† †Indian Institute of Technology Madras, Chennai-600036, India. ‡Current address: University of Alberta, Edmonton, Alberta T6G1H9, Canada E-mail: [email protected] Phone: +91-44-2257-4177 Abstract Input design is the process of designing informative inputs which maximizes the information about system dynamics using minimal cost and resources in a system identification experiment. In this work, we focus on designing inputs based on the economic framework for identification of constrained processes. The process is nominally constrained and feasibility during identification is ensured by “backing-off” (or moving away) from the active constraints. The economic cost or penalty due to the “back-off” is minimized while ensuring feasibility and satisfying constraints on the quality of the identified model. The proposed economic based formulation results in a non-convex optimization problem. This problem can be solved using bisection algorithm for SISO systems. We present a two-stage iterative solution algorithm to solve for MIMO systems where a convex optimization problem is solved at each stage. The optimal input is realized as a white noise passing through an M-tap filter. The filters coefficients are obtained by a non iterative spectral factorization technique. The efficacy of the proposed approach has been successfully demonstrated in a mass spring damper system as well as in a Van der Vusse reactor system.

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

Keywords:

System identification, economical input design, input design, convex

optimization, robust constrained, semidefinite programming.

Introduction Optimal input design for system identification has been an intensive research area for the last four decades. The objective of any identification exercise and hence, input design is to estimate parameters in a dynamic model efficiently and precisely. The choice of input or an “optimal input design" is the most fundamental step in identification in the sense that it may determine the nature and accuracy of the system characteristic from the identification procedure. 1 Typically, the input design problem is posed as an optimization problem where the cost function is some scalar measure of the Fisher information matrix subject to constraints on inputs and outputs. Based on the scalar measure of cost function, input design problem can be classified as A, E, or D-optimal design problem. 1–3 The input design can be addressed in the Bayesian framework, 4–6 or robust framework. 7 In the Bayesian approach, the problem is formulated in a probabilistic framework, where the objective is to maximize expected utility or minimize risk. The utility function is chosen based on the goal of the experiment design. The robust approach formulates a min-max problem, where it minimizes the worst case scenario over an uncertain set of parameters. 7 The renewed interest in input design was spurred by developments in convex optimization theory and algorithms, in particular, semidefinite programming. 8 Some of the classical and recent input design problems can be formulated as convex optimization problems. In particular, it is often posed as Semidefinite programming problems (SDPs) or generalized eigenvalue problems (GEVPs). These problems can be solved using available convex optimization tools such as CVX 9 or YALMIP. 10 However, converting input design problem into SDPs are not often trivial. Recently, the least costly identification framework has been introduced where the weighted input and output power is minimized subject to constraints on model quality which is represented as a function

2


Page 2 of 41

Page 3 of 41 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


of parameter covariance matrix. 11 This can be interepreted as a “dual" to the conventional input design problem. Recently, there has been increased research interests in designing optimal input for non-linear system. A robust approach of input design for non-linear system is presented using graph theory approach where a scalar cost function of information matrix is minimized over a set of marginal distributions of stationary processes. 12 A similar notion of cost function is also derived to design an optimal input for non-linear FIR systems. 13 In practice, the identification of process plants and equipments is carried out on a running plant and any perturbation in input will introduce additional variability in the output. 14 Hence, designer has to respect certain constraints on inputs and outputs to ensure safety, feasibility and product quality during the identification experiment. This is popularly known as “plant friendly identification”. In process industries, economics often stipulate that the plant or equipment be operated at or near constraints to maximize profit. 15,16 However, due to modeling error, parametric uncertainty, unexpected disturbance or additional variability introduced in the input and output due to the identification experiment, this will result in higher risk of constraint violation. Hence, the experimental conditions such as the operating point and the perturbation signal must be chosen so that the risk of constraint violation is within satisfactory limits. One solution is to “back-off" from the constraints (or equivalently moving away from the constraints inside the feasible region) which results in a trade-off between profitable operation and risk of constraint violation. Simultaneously, the overall goal of the experiment (in terms of quality of the identified model) has to be satisfied. This is similar to the least costly framework, except that we quantify the cost of the experiment in terms of the deviation from the nominal operating policy (or equivalently, deviation from constraints) rather than input or output variance. We address this problem by first describing the expected dynamic operating region (EDOR) of the plant described by LTI models. 17,18 The EDOR is the region in the input/output or state space where the system is likely to operate (modulo a probability level). 16 This is a function of the disturbances and additional perturbation introduced due to iden-

3



Page 4 of 41

tification and under fairly general conditions, it is an ellipsoid. The aim is to design an appropriate identification experiment so that the cost of the experiment (quantified in terms of deviation from optimality) is minimized. This is carried out subject to the overall goal of the experiment being satisfied and the dynamic operating region being feasible. The decision variables are the experiment conditions, viz., the operation point around which the perturbation has to be carried out and the input spectrum. The input spectrum is finitely parameterized by requiring it to be realizable by a M-tap FIR filter. The filter coefficients are obtained by a non-iterative approach for spectral factorization. 19 A two step iterative method has been proposed to solve the resulting non-convex problem. 17,18 In this work, we present a bisection solution algorithm for designing inputs in SISO systems which is guaranteed to converge to the globally optimal solution. In addition, we explore the influence of filter length on the back-off value and EDOR. This paper is organized as follows: Firstly, we present the preliminaries of system identification followed by economic based framework for input design problem. Next, the optimization formulation for input design is presented and then computationally efficient solution algorithms are proposed. Finally, illustrations are provided to demonstrate the proposed input design approach.

PRELIMINARIES In this section, we briefly present the preliminary aspects of system identification. Let us consider the identification of a discrete-time, stable, linear time invariant system operating in open loop. The system is modelled as

M : y(k) = G(q −1 , θ)u(k) + H(q −1 , θ)e(k)

(1)

where y ∈ Rny and u ∈ Rnu are output and input signals respectively. e ∈ Rny denotes an additive Gaussian white noise sequence with zero mean and variance Λ = diag[λ1 , λ2 · · · λny ]. 4


Page 5 of 41 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


q −1 is a back-shift operator (that is, q −1 u(k) = u(k − 1)) and, G(q −1 , θ) and H(q −1 , θ) are stable rational transfer function matrices of dimensions ny × nu and ny × ny respectively. We assume that the identified model is flexible enough to capture dynamics of the true system S, i.e., the true system is in the model set. 20

S : y(k) = G0 (q −1 )u(k) + H0 (q −1 )e(k)

(2)

In addition, we make the following assumptions: • There are no pole-zero cancellations. • The input is a stationary signal with power spectrum Φu (ω) ∈ Rnu ×nu defined on [−π, π]. R[τ ] ∈ Rnu ×nu is the matrix valued correlation coefficient of input u at lag τ and forms a Fourier transform pair with joint input spectrum Φu (ω). • H(0, θ) = Iny The model parameters θ are identified within the prediction error framework where the objective is to minimize one step ahead prediction error (k, θ), 21

θˆ = arg min θ

T −1 (k, θ) Λ (k, θ) k=1

P N 1 N

(k, θ) = y(k) − yˆ(k|k − 1, θ)

(3) (4)

where N is the number of data points, yˆ(k|k − 1, θ) is the one step ahead prediction of y which is computed from the assumed model. Under some mild conditions, the estimated parameter follows a normal distribution: 20 √ N →∞ N (θˆ − θ0 ) −−−→ N (0, Pθ )

(5)

where Pθ is a covariance matrix which characterizes the uncertainty in the parameter esti-

5



Page 6 of 41

mates, and it can be expressed as, "

∂(k, θ) −1 ∂(k, θ) T Pθ = E Λ ∂θ ∂θ

#−1 (6)

For an asymptotic unbiased estimator, the covariance of the parameter estimate is given by the expression Pθ = Mθ−1

(7)

where Mθ is the Fisher information matrix which can be further partitioned as input dependent matrix Mβ and input independent matrix Mγ as follows: 

  Mθ = 

Mβ (u, θ)

Mue (u, Λ, θ)   T (u, Λ, θ) Mγ (Λ, θ) Mue

(8)

For an open loop system, Mue = 0, because noise and inputs are uncorrelated with each other. The individual entry (j, k)th of the partitioned matrix Mβ and Mγ can be calculated as follows:. 22 1 = 2π

[Mβ ]j,k

Zπ

˙ ˙∗ tr{Φ−1 v Gk Φu Gj }dω

(9)

−π

[Mγ ]j,k

1 = 2π

Zπ

˙ ˙∗ tr{Φ−1 v Hk ΛHj }dω

(10)

−π

where G˙ j , H˙ j are the partial derivatives of the G and H matrices respectively with respect to the j th parameter. Superscripts (*) present on G and H denote conjugate transpose of the particular matrix. Φu is the joint input spectrum and Φv = HΛH ∗ is the noise spectral matrix. Thus, Mθ can be partitioned into a block diagonal matrix such that the upper block (Mβ , u) is influenced by the input and the lower block is not. Hence, in input design problems, it is pertinent to consider only the input dependent part Mβ .

6


Page 7 of 41 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


PROBLEM FORMULATION In this section, we will briefly present the traditional input design problem followed by the economical input design formulation.

Optimal input design In a traditional input design formulation, the objective is to minimize some scalar norm of information matrix subject to constraints on signal power, product quality or combination of both. 2,3,23 The general form of formulation is presented as follows:

min Φu

s.t

J (Φu )

signal power constraint,

(11)

quality constraint, or combination of both

where Φu ∈ Rnu ×nu is an input spectrum. The diagonal elements of Φu signify the spectrum of j th input (Φuj ) and the off-diagonal elements (Φujk = Φukj ) denote cross spectrum between j th and k th inputs. The most common objective functions (J ) are: − log det (Mβ ) (Doptimal), −λmax (Mβ ) (E-optimal), tr(Mβ−1 ) (A-optimal), tr(W Mβ−1 ) (L-optimal) etc. In the least costly framework, the weighted sum of the input and output power is minimized while ensuring that the quality of an identified model is sufficiently accurate for the intended use. 11 For example, the E-optimal formulation can be modified in the least costly framework as follows:

min αu tr(Σu (ω)) + αy tr(Σy (ω)) Φu

s.t. Mβ (ω) ≥ ηI

7


(12) (13)


Σu (ω) and Σy (ω) denote the respective input and output power and these are functions of input spectrum Φu (ω). αu and αy are respective weights on the input and output power. ηI is an upper bound on the covariance matrix of the estimated parameters and is a user defined parameter which controls the quality of the identified model. A higher value of η implies that the identified model have low parameter variance and will require more identification effort in terms of input excitation. More general control relevant constraints can also be imposed. 11,24 The choice of weights is based on the relative importance or costs associated with the respective variations. The implementation depends on whether there are constraints or not. E.g., if the system is unconstrained, the optimum input can be safely generated and the system perturbed with the input. However, if there are constraints, On the other hand, if the nominal operating point is at the limiting values of inequality constraints, then the designed input might result in constraint violations due to the presence of measurement errors and disturbances, etc. Hence, there is no guarantee that the above approach will result in a feasible operation during the experiment. Further, it is not clear how to choose the weights appropriately. Therefore, the focus of this work is an input design formulation for constrained systems where the above concerns are addressed in an economic framework.

Economical input design Economic input design for identification deals with the problem of determining the profitable operating point and optimal input signal while ensuring feasibility during identification and maximizing the quality of the identified model. 17 Given an economic cost function that is a function of the inputs and outputs, J(u, y), the optimal operating point [u∗T , y ∗T ]T of the plant is determined by solving a static, non-linear optimization problem by minimizing J(u, y) subject to the plant model and constraints. These constraints are usually in the form of input, output or resource constraints. Often the optimal operating point (OOP), [u∗T , y ∗T ]T is at the intersection of constraints, i.e., some constraints are active as shown in Figure 1. Uncertainties in the form of disturbances can cause constraint violations. In addi8


Page 8 of 41

Page 9 of 41 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


𝑦 OOP

𝑢

𝑦𝑚𝑎𝑥

BOP

𝑢𝑚𝑖𝑛

𝑢𝑚𝑎𝑥

𝑦𝑚𝑖𝑛

Figure 1: Basic idea of profitable input design at the constrained optimal point in the presence of uncertainties tion, if we intend to carry out system identification experiments at this profitable operating point, the risk of constraint violation increases because of input and output perturbations. Typically, an informative experiment requires sufficient excitation. Hence, a conservative solution is to “back-off” or equivalently, operate the plant farther inside the feasible region T T ] , which will result in non-profitable operat the backed-off operating point (BOP), [uTss , yss

ation. 16,25,26 Therefore, it is necessary to trade-off the profitability and the objectives of the experiment. Hence, in such situations, a more relevant cost associated with the identification experiT T ment is the loss in operating profit incurred by operating the plant at [uTss , yss ] rather than

at [u∗T , y ∗T ]T . Therefore, the aim of this work is to determine the experiment conditions that minimizes this loss. These include the operating point around which the identification experiments is carried out and the input spectrum. This is formulated as optimization problem where the objective is to minimize the loss in operating profit such that the input design can be carried out without violating the constraints. In other words, the dynamic operation is feasible during the identification exercise. Define the variables, z = [uT , y T ]T ∈ Rnz , nz = nu + ny . Let us rewrite all the constraints 9



Page 10 of 41

in the form of aTj z ≤ bj , j = 1, . . . , nc . Also, we assume a linear form of the cost function J(uss , yss ) which is a good approximation at the constrained optimal point. Since the disturbances are random, we consider a stochastic framework. 16–18 We define the expected dynamic operating region (EDOR), E ∈ Rnz as the region in the z space such that the probability of z ∈ E is sufficiently high. 16,26 . Now, the optimization formulation can be expressed as

min

Φu ,uss ,yss

JuT (uss − u∗ ) + JyT (yss − y ∗ )

(14)

s.t. yss − y ∗ = K(uss − u∗ )

(15)

Mβ (Φu (ω)) ≥ ηI

(16)

aTj z ≤ bj ∀z ∈ E, j = 1, . . . , nc

(17)

Φu (ω) 0, ω ∈ [−π, π]

(18)

In the above formulation, we quantify the economic cost of identification in terms of deviation from the optimal operating point. The constraints are operational constraints (operating point should lie on steady state operating line), model quality constraints, and, input and output constraints. The objective function and the constraints are linear in terms of decision variables. However, the problem is hard to solve because of the following reasons: • Since Φu is defined over a continuous frequency ω ∈ [−π, π], the corresponding constraint (18) is infinite dimensional. • The constraint (17) involves a finite number of variables, but infinite number of con straints, since aTj z ≤ bj should be satisfied for all z ∈ E. This enforces the expected dynamic region (corresponding to the specified probability level) to be within the constrained operating space defined by aTj z ≤ bj .

10


Page 11 of 41 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


Finite dimensional formulation In this section, we reformulate the problem and show how it can be converted into a finite dimensional problem. First, we discuss the finite parameterization of input spectrum. Next, we characterize the EDOR and present a finite dimenstional equivalent of (17).

Parameterization of input spectrum The information matrix depends on the input spectrum, which is infinite dimensional. Therefore, in this section, we present the finite parameterization of input spectrum to deal with semi-infinite constraint (18) presented in the economic input design formulation. A common method of finitely parameterizing the specturm of the j th input in terms of a finite number (m) of basis functions as follows: 23 (M −1)

Φj (ω) =

X

cj (τ )e−iτ ω ≥ 0 ∀ ω ∈ [−π, π]

(19)

τ =−(M −1)

where i =

√

−1. The cross-spectrum between j th and k th input is given by 27 : (M −1)

Φjk (ω) =

X

cjk (τ )e−iτ ω

(20)

τ =−(M −1)

Using the above parameterization the joint input spectrum is written as: 21 (M −1)

Φu (ω) =

X

C(τ )e−iτ ω 0 ∀ ω ∈ [−π, π]

τ =−(M −1)



 c1 (τ ) c12 (τ )   c12 (τ ) c2 (τ )  C(τ ) =   ··· ···   c1nu (τ ) c2nu (τ )

11

 c1nu (τ )   · · · c2nu (τ )    ··· ···    · · · cnu (τ ) ···


(21)


Page 12 of 41

where C(τ ) is a coefficient matrix of dimension Rnu ×nu and must be such that C(τ ) = C(−τ ). 21 Instead of parameterizing the full spectrum, we may equivalently work with a parameterization in terms of Ψ(eiω ) as: 23,28

Φu (ω) = Ψ(eiω ) + Ψ∗ (eiω ) 0, ∀ ω ∈ [−π, π] M −1 C(0) X lω + C(τ )e−iτ ω Ψ(e ) = 2 τ =1

(22) (23)

We can interpret the above parameterization as restricting the input signal to be white noise sequences passing through an FIR filter of length m, and hence C(τ ) can be interpreted as the correlation sequence (R(τ ) = E u(k)u(k − τ )T ∈ Rnu ×nu ) of the input signal. 21 In this case, the input spectrum (Φu ) and Ψ( eiω ) can be expressed as: 23,28

Φu (ω) = R(0) + 2

m−1 X

R(τ )cos (τ ω) 0 ∀ ω ∈ [−π, π]

(24)

τ =1 m−1

R(0) X Ψ(e ) = + R(τ )e−iωτ 2 τ =1 iω

(25)

Since Ψ is complex, it can be interpreted as a transfer function of a finite dimensional system. In order to obtain a state space realization of Ψ, we need to define Aφ , Bφ , Cφ , Dφ in terms of R(τ ), τ = 0, . . . , m − 1 such that

Cφ (eiω I − Aφ )−1 Bφ + Dφ =

R(0) + R(1)e−iω + · · · R(M − 1)e−iω(M −1) 2

12


(26)

Page 13 of 41 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


One possible choice of the above representation is controllable canonical form:       Aφ =      

0

0

0 0

···

0

···

Inu

0

0 .. .

Inu 0

0 . 0 · · · .. . . . . .. . . .

0

0

0

Inu 0





          

     Bφ =      

0 0 .. . 0

(M −1)nu ×(M −1)nu

Cφ = [R(1), R(2), . . . , R(m − 1)]nu ×(M −1)nu

Inu

Dφ =

            (M −1)nu ×nu

1 [R(0)]nu ×nu 2

It can be easily verified that Ψ(ω) = (Cφ (eiω I − Aφ )−1 Bφ + Dφ ). This realization is not unique and a simple matrix transformation realizes the same transfer function. 28 For example, consider a case where nu = 2 and M = 3, i.e., the input spectrum is finitely parameterized using three parameters R(0), R(1), R(2), where R(τ ) ∈ R2×2 , τ = 0, 1, 2. For this choice of parameterization, we have:

Ψ(eiω ) =

R(0) + R(1)e−iω + R(2)e−ilω 2

(27)

A state space realization is as follows:



 0 0 0   0 0 0  Aφ =   1 0 0   0 1 0





0   1 0    0 1 0     Bφ =    0 0 0     0 0 0

        

1 Cφ = [R(1) R(2)] ; Dφ = R(0) 2 The positivity condition of input spectrum (Φu 0) can be imposed by invoking the following Kalman-Yakubovich-Popov (KYP) lemma. 23,28

13



Page 14 of 41

Lemma 1. Given a discrete-time linear system (Aφ , Bφ , Cφ , Dφ ), Aφ stable, (Aφ , Bφ , Cφ ) minimal and Dφ + DφT ≥ 0. The transfer function Ψ(ω) = Cφ (eiω I − Aφ )−1 Bφ + Dφ satisfies Φu (ω) = Ψ(eiω ) + Ψ∗ (eiω ) 0 ∀ ω ∈ [−π, π]

if and only if there exists real symmetric matrix Qφ such that the following linear matrix inequality is satisfied 

 ATφ Qφ Aφ

CφT

ATφ Qφ Bφ

−  Qφ −  0  T T T Cφ − Bφ Qφ Aφ Dφ + Dφ − Bφ Qφ Bφ

(28)

where Qφ ∈ S(m−1)nu ×(m−1)nu . For proof, the reader is referred to. 23,28,29

Characterizing the expected dynamic operating region In this section, we characterize the EDOR and transform the infinite set of constraints (17) into an equivalent representation of finite set of constraints using convex optimization techniques. Given the first and second moments of z in terms of as follows,

Ez = zss

(29)

E(z − zss )(z − z˜ss )T = Σz (ω)    Σu (ω) Σuy (ω)  Σz (ω) =   Σyu (ω) Σy (ω)

(30)

14


(31)

Page 15 of 41 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


where Σu , Σy and Σyu denote the input variance, output variance and covariance respectively. We then express them in the frequency domain as follows:

Σu

1 = 2π

Σy =

1 2π

Zπ −π Zπ

Φu (ω)dω

(32)

(GΦu G∗ + HΛH ∗ )dω

(33)

−π

Σyu = Σuy

1 = 2π

Zπ

G(eiω , θ)Φu (ω)dω

(34)

−π

Recall that the noise ek is assumed to be Gaussian and the input is obtained by filtering a Gaussian white noise signal through an m-tap FIR filter. Therefore, the vector z is jointly Gaussian and is completely characterized by the first and second order moments. Hence, the EDOR is an ellipsoid and for a specified probability level, it can be completely characterized in terms of the first and second order moments of z as follows:

E = {z|z = zss + αΣ1/2 ˜, ||˜ z || ≤ 1} z z

(35)

where α depends on the probability level. Theorem The set of linear constraints aTj z ≤ bj , z ∈ E can be expressed as a set of second 1/2

order cone constraints of the form aT zss + αkΣz aj k2 ≤ bj . 1/2

Proof: 30 Recall that z := zss +αΣz z˜, ||˜ z || ≤ 1. Consider the infinite dimensional constraint, aTj z ≤ bj ∀z ∈ E

(36)

⇐⇒ sup{aTj z | z ∈ E} ≤ bj

(37)

⇐⇒ aTj zss + sup{aTj αΣ1/2 z z | kzk2 ≤ 1} ≤ bj

(38)

⇐⇒ aj T zss + αkΣ1/2 z aj k2 ≤ bj

(39)

Remark 15



Page 16 of 41

The above second order cone constraint can also be represented as an LMI using S-procedure,  T  −κj − aj z˜ss + bj  1/2 ( α2 aT Σz )T

α T 1/2 a Σz 2

κj I

   0; κj > 0

(40)

The constraint (35) can also be interpreted as robust version of uncertain constraints aTj z ≤ bj , where z is a random variable with mean zss and covariance Σz . Now, the finite dimensional formulation of the economically optimal input design problem can be presented as follows:

min R(0),...,R(m−1),uss ,yss ,Q

JuT (uss − u∗) + JyT (yss − y ∗ )

s.t.yss − y ∗ = K(uss − u∗ )    Σu (R) Σuy (R)  Σz (R) :=   Σyu (R) Σy (R) Mβ (R) ≥ ηI T α||Σ1/2 z aj || + aj zss ≤ bj , j = 1, . . . , nc   T CφT − ATφ QBφ   Q − Aφ Q φ Aφ 0  T T T Cφ − Bφ Qφ Aφ Dφ + Dφ − Bφ QBφ

(41)

Q = QT

where the last LMI (41) follows from the KYP lemma and input parameterization described in Section 3.1. The above formulation is finite dimensional and the decision variables are input autocorrelation sequences and the backed-off operating points. However, the presence of Σ1/2 in the constraint makes the problem non-convex. To address this issue, we propose two solution methodologies.

16


Page 17 of 41 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


Solution Techniques In the previous section, we formulated the input design problem that will result in a profitable and feasible operation while guaranteeing that the identified model is sufficiently precise. For this purpose, we characterized the dynamic operation of the input design experiments using ellipsoids. The center of the ellipsoid defines the profitable operating point around which the system is perturbed. However, the formulated plant friendly input design problem is non-convex and hence the above formulation is computationally intractable due to the presence of nonlinear matrix equality constraints. Equivalently, the problem can be interpreted as determining the center zss and the orientation (decided by the excitation) while ensuring feasibility of the ellipse and model quality. Therefore, the focus of this section is to present computationally efficient solution techniques. 17,26 First, we describe a novel two-stage iterative algorithm that is applicable for Multi-Input Multi-Output systems. Next, we present a bisection algorithm that can be used to determine globally optimal solutions for Single-Input Single-Output systems.

Two-stage algorithm As mentioned previously, the optimization formulation aims at determining the nearest backed-off operating point to the optimal point such that the input design experiments are feasible. In other words, the algorithm should find the economic back-off point that has the least variability in actively constrained variables. To this end, we follow a novel two-stage iterative solution algorithm. 17,26 The basic idea in the first stage is to find the smallest (in terms of trace) feasible ellipsoid Σz . For ease of exposition, we assume that the constraints are in the form of bound constraints, i.e., zj ∈ [dmin (j), dmax (j)]. Other arbitrary linear constraints of the form aTj z ≤ bj can be addressed with appropriate modifications. To

17



Page 18 of 41

ensure feasibility, we impose the following constraints on the individual variances,

Σz (j, j) = σz2j ≤

1 (dmax (j) − dmin (j))2 j = 1, · · · (nz ) 4α2

(42)

where σz2j is the variance of zj It is important to note that, in a traditional input design problem, the upper limit is assumed to be provided by the designer whereas, in the economical framework, it is defined in terms of the available feasible space given by the variable bounds (dmin , dmax ). To this end, we define the following constraints with respect to variance of the 2 , k th variable, σz,k

2 σz,j ≥

2 δjk σ2 α2 z,k

(43)

where the iterative parameters δjk are chosen such that the BOP selected in stage 2 is used to select the new minimum variance ellipsoid that forces the BOP close to OOP. The parameter δjk is defined as

δjk =

distance of variable j f rom its closest bound distance of variable k f rom its closest bound

(44)

Initially the parameter δjk is set as zero, indicating the lower bound on the variances to be positive, which is physically meaningful. The Stage 1 optimization problem is given below

18


Page 19 of 41 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


min

tr(Σz )

s.t.

Mβ (R) ≥ ηI   T CφT − ATφ QBφ   Q − Aφ QAφ 0  T T T Cφ − Bφ QAφ Dφ + Dφ − Bφ QBφ   Σuy (R)   R[0] Σz =   0; Σyu (R) Σy (R)

R,Q

Σz (j, j) ≤ 2 ≥ σz,j

1 4α2

2 δjk σ2 , α2 z,k

(dmax (j) − dmin (j))2 , j = 1, . . . , nz j 6= k

Q = QT Solution of Stage 1 results in a feasible covariance ellipsoid Σz . This is used to find the approximation to the backed-off operating point in stage 2 as follows. In the second stage, this covariance ellipsoid is used to determine the closest possible backed-off operating point (BOP), zss , to the constrained optimal operating point OOP, z ∗ such that the input we design will be feasible. The Stage 2 optimization problem can be formulated as min

Ju T (uss − u∗ ) + Jy T (yss − y ∗ )

s.t.

yss − y ∗ = K(uss − u∗ )

uss ,yss

1/2

α||Σz aj || + aTj zss ≤ bj , j = 1, . . . , nc The backed-off operating point obtained after the first iteration might be sub-optimal. For nominally constrained processes, the surface of the optimal dynamic operating region should touch the boundary of the feasibility window. Therefore, the current information from the second stage (i.e., BOP) is used to create lower bounds on the variances by defining the parameter δ. This parameter defines the closeness of the BOP obtained at the current iteration to the OOP and this information is used for bounding the individual variances to recompute Σz in the first stage. At the start of the procedure, δ is initialized to 0. It is 1/2

important to note that Σz

is not a decision variable since Σz is known from first stage.

Now it can be easily recognized that both the stages contains only convex constraints which could be easily solved using CVX, a package for specifying and solving convex programs. 9 19



iter−1 iter Given two successive iterates, zss and zss this process is iterated until the convergence iter-1 iter k2 ≤ is satisfied where being the prescribed tolerance limit. The − zss criteria kzss

detailed solution algorithm is presented in Table 1. The proposed solution methodology is quite general and applicable for both SISO and MIMO systems. However, the theoretical properties of the algorithm, viz., convergence to the global optimum are still under investigation. However, an alternate solution approach using bisection algorithm is proposed for SISO systems which is guaranteed to converge to the global optimum. Table 1: Two-stage algorithm to determine the economic back-off operating point for input design 1 Initialize the parameter δ = 0. 2 Find Σz by solving the following convex problem (Stage 1), min tr(Σz ) R,Q

s.t.

Mβ (R) ≥ ηI Q − ATφ QAφ CφT − ATφ QBφ 0 Cφ − BφT QAφ Dφ + DφT − BφT QBφ Q = QT R[0] Σuy (R) Σz = 0; Σyu (R) Σy (R) Σz (j, j) ≤ 4α1 2 (dmax (j) − dmin (j))2 ; j=1,· · · nz 2 δjk σ2 , α2 z,k

2 ≥ σz,j

j 6= k

If no feasible Z can be found, exit. 1/2 3 Compute P = Σz . Find the BOP (zss ) by solving the following convex problem (Stage 2), min

Ju T (uss − u∗ ) + Jy T (yss − y ∗ )

s.t.

yss − y ∗ = K(uss − u∗ ) α||P aj || + aTj zss ≤ bj

uss ,yss

4 Terminate on convergence. Otherwise, update δ using (44) from the current value of BOP and proceed to Step 2.

20


Page 20 of 41

Page 21 of 41 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


𝑦 OOP

𝑦𝑚𝑎𝑥

𝑢

𝑦𝑏𝑖𝑠 𝑦𝑏𝑖𝑠 =

𝑦𝑂𝑂𝑃 + 𝑦𝐵𝑂𝑃 2

BOP

𝑦𝐵𝑂𝑃 𝑢𝑚𝑎𝑥

𝑢𝑚𝑖𝑛

𝑦𝑚𝑖𝑛

Figure 2: Illustration of Bisection method

Bisection algorithm In this section, we present an alternative solution technique for solving the proposed input design problem. As noted previously, the operating point and size and orientation of the covariance ellipsoid play a vital role for profitable input design. For single input single output systems, the backed-off operating point can be uniquely determined by uss ∈ R. The basic idea of the bisection algorithm is to perform a feasibility check of the covariance ellipsoid for a guess uss ∈ R. If feasible, we move towards the BOP along the steady state line, else in the opposite direction. For this purpose, we set δ = 0 and replace the individual variance constraints (42) of the Stage 1 problem with the following constraints

σu2

Economically Optimal Input Design Approach for Identification of

Recommend Documents