Identification of MIMO Continuous-Time Models Using Simultaneous

Article pubs.acs.org/IECR

Identification of MIMO Continuous-Time Models Using Simultaneous Step Inputs Salim Ahmed* and Syed A. Imtiaz Department of Process Engineering, Memorial University, St. John’s, Newfoundland and Labrador A1B 3X9, Canada ABSTRACT: A new approach for identification of multi-input multi-output (MIMO) continuous-time transfer function models using simultaneous step input signals is proposed. MIMO processes exhibit directionality which implies that the output gains depend on the input magnitudes as well as the ratios of the inputs. Due to this directionality issue, MIMO models estimated by sequentially changing one input at a time often do not result in satisfactory tracking performance when used for model based control. In the proposed methodology all of the inputs are changed simultaneously to resemble controlled conditions. While the identification tests remain multi-input in its true sense, the parameter estimation steps involve estimation of the parameters of a single transfer function at a time. Moreover, the time delays are estimated in the same way as the coefficients in the model. Simulation results are presented to demonstrate the robustness of the methodology and its applicability to processes with high input−output dimensions. Identification and control results of a simulated distillation column are also presented; a dynamic matrix controller (DMC) results in better control performance with the model estimated using the proposed methodology than with the model estimated using sequential inputs.

1. INTRODUCTION 1.1. Background and Motivations. A unique characteristic of multi-input multi-output (MIMO) systems is their directionality, meaning that the magnitudes of the outputs are not only dependent on the magnitudes of the inputs but also on their ratios (often called input directions). For illconditioned systems, directionality poses challenges in both identification and control of the system, especially in the weak directions. In order to perform effective control in the weak directions, the corresponding models are required to be very precise. Plant−model mismatch resulting from inadequate excitation in the weak directions can lead to significantly poor control performance.1,2 Our motivation for working on this problem comes from the challenges in model predictive control (MPC), especially in the linear programming (LP) layer where the gains of the process models are used to calculate the control moves. A binary distillation column is used as an example to illustrate the control problem. This system has been extensively studied in the literature.3,4 The nonlinear aspects of distillation columns and the associated difficulties in their identification have also been highlighted by Luyben5 and Martin et al.6 The control objective for distillation is to maintain nearly pure lighter component as distillate and nearly pure heavier component in the bottoms product. There can be many different control configurations for the system depending upon the types of controllers in use. The two controlled variables (CVs) for this system are the distillate concentration and the bottoms concentration. A common control configuration is to use the reflux rate and the reboiler duty as the manipulated variables (MVs). A vector plot depicting the relative directions of the inputs and outputs is shown in Figure 1. As shown in the figure, the effects of reboiler duty and reflux on the CVs are almost opposite to each other. Although the control problem is difficult with the multivariate configuration, it is often selected © XXXX American Chemical Society

Figure 1. Vector plot representing the MVs for a desired change in CVs.

as it gives the controller more degrees of freedom. Typically in binary distillation, the control objective is to recover both top and bottom products in their pure forms. Such a control objective is nearly orthogonal to both MVs. Thus, this control direction is considered a weak direction relative to the MVs. In order to gain a small improvement in both concentrations, large changes have to be made in both of the MVs. A numerical example will further highlight the issues involving control and identification. Consider a distillation column with the model gains given in Table 1. For this column, if a change of 0.1 mol % is desired in both top and bottom Received: February 3, 2015 Revised: June 24, 2015 Accepted: July 2, 2015

A

DOI: 10.1021/acs.iecr.5b00481 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

constructing the frequency response using the fast Fourier transform (FFT) and inverse FFT. Beside the state space approach, MIMO systems are also identified as nonparametric frequency responses as in ref 16. Another group of methods17−19 first obtains the frequency response and subsequently estimates the step response of the individual input−output channels and the transfer functions. For the purpose of simulation, prediction, or controller design, the transfer function models are more accurate and need shorter test time than the nonparametric frequency response models.20 The procedure to obtain the transfer function from a frequency response is not straightforward. Jin et al.21 proposed a closed loop identification method for multivariable continuous-time models using a nonlinear stochastic direct search optimization algorithm. Other work using optimization methods are refs 22 and 23. The closed loop identification method for multivariable processes from step responses has also been proposed using the B-Spline expansion.24 Wang et al.25 used the relay feedback method for identification of MIMO systems. An iterative nonlinear separable least-squares method that requires discretization of the continuous-time model was proposed in ref 26 for multi-input systems. Most of these methods use sequential step test data, where step inputs are applied to manipulated variables one at a time while keeping the other inputs as constants; this process is then repeated for each of the inputs. For example, ref 27 discussed open loop identification of MIMO continuous-time models using rectangular pulse responses. A genetic algorithm based closed loop identification method was proposed in ref 28. A sequential loop closing approach has been used for identification of MIMO models for autotuning29 and for iterative identification of ill-conditioned processes.30 An algorithm to identify a two-by-two transfer function model from closed loop sequential step testing has been proposed in ref 31. The importance of MIMO identification has been highlighted in ref 32. Several authors have proposed the idea of simultaneous excitation.20 The importance of MIMO input design and a comparative analysis of the performance of the estimated model used for model predictive control has been presented in ref 33. Issues with closed loop identification of MIMO systems have been discussed in ref 34. Use of the optimization method for parameter and input estimation in the case of incomplete output data was discussed in ref 35. A modified version of the bee colony algorithm for structural system identification is presented in ref 36. The above review highlights the significance as well as the lack of development in continuous-time identification for direct estimation of MIMO model parameters from truly multi-input test data. In this work we propose a new methodology to directly estimate continuous-time transfer function model parameters from simultaneous step test data. While the identification tests remain multi-input in its true sense, the parameter estimation steps involve estimation of the parameters of a single transfer function at a time. Moreover, the time delays are estimated in the same way as the coefficients in the model. The contribution of this work is tied to the above unique aspects of the proposed approach. The remainder of the article is organized as follows. The identification methodology is detailed in Section 2. Identification results of a simulated distillation column is presented in Section 3 followed by simulation results of higher dimension models in Section 4. Concluding remarks are presented in Section 5.

Table 1. Steady State Gains between CVs and MVs for a Distillation Column CV/MV

bottoms impurity (mol %)

distillate impurity (mol %)

reflux rate (kbpd) steam rate (Mlb/h)

0.46 −0.61

−0.18 0.25

products, and this will require a change of 3.7 kbpd in the reflux flow and 2.9 Mlb/h in the reboiler duty. Now assume a 10% model error in the gain between reflux and bottom composition due to the inaccuracies in the identification process. Using the new gains the required MV moves are 7.6 kbpd and 5.5 Mlb/h for the same desired changes in the output. This is twice the MV changes of that actually required to achieve the quality target changes. This error is very large considering the small changes desired for the quality target. Therefore, appropriate methods are required for step testing of ill-conditioned systems to correctly identify these models. The traditional step testing approach of moving one MV at a time, e.g., in the case of a distillation column, moving only reflux or only reboil, will change the cut point in the column (the so-called strong direction). But this will make little or no effect on separation or the mass transfer between the downcoming liquid stream and the vapor stream going toward the top of the tower (the weak direction). Unless both MVs are moved simultaneously, the mass transfer/separation mechanism will not be in effect and the model gain will be underestimated. In a control scenario both MVs will be moved simultaneously to get pure components in top and bottom. Therefore, the identification step should mimic similar conditions so that the identified model captures the effect of both mechanisms. Of course, this raises the valid issue of how MV cross-correlations will affect the results. It is well-known that if the MVs are highly correlated, the model identification results can be adversely affected. Hence identification methods need to be devised that can handle the correlation between the input variables. 1.2. Brief Review of the Literature. Identification of an ill-conditioned system requires that the inputs are designed in such a way that the weak direction dynamics of the system are fully excited. Sequential steps applied to different input channels cannot excite the process in the weak direction. In this scenario simultaneous changes (i.e., correlated inputs) of the inputs can give the most relevant model for control. In discrete time system identification this problem has been dealt.7 The singular value decomposition (SVD) based methods use correlated inputs with input directions aligned with the eigenvectors. It allows a maximum gain in the output direction. However, there is not much work related to this issue in the field of continuous-time identification, although there are advantages of continuous-time MIMO model identification.8,9 The theoretical development of continuous-time MIMO system identification has mainly been based on the state-space framework.10−13 Aziz et al.14 proposed a two-stage approach for identification of MIMO continuous-time models that first estimates the step response using subspace approaches; the step response is then used to estimate the model parameters. Jeng and Lin15 proposed closed loop identification techniques for multivariable processes to simultaneously estimate the process and disturbance models by first estimating the finite impulse response model using a subspace algorithm and then by B


Article

Industrial & Engineering Chemistry Research Nu

2. PROPOSED METHODOLOGY The basic idea behind the proposed approach is to decompose the MIMO data matrix into equivalent single-input singleoutput (SISO) data matrices by performing a series of algebraic manipulations. The main advantage of this decomposition comes from the fact that, in the parameter estimation steps, parameters of each transfer function are estimated separately by solving independent equations. Thus, the complexities in the parameter estimation step for MIMO systems, that arise due to the large numbers of parameters, are avoided. The mathematical formulation for the proposed methodology is described below. Consider the MIMO system in open loop described in Figure 2 where uj with j = 1, 2, ..., Nu are the inputs; xi with i =

Xi , l =

∑ Gi ,jUj ,l

with

i = 1, 2, ..., Ny (3)

j=1

Note that the second subscripts in Xj,l and Uj,l refer to the lth number of experiment. Equation 3 can be written for the ith output and for all the experiments, i.e., for l = 1, 2, ..., Nu, and can be combined and represented in a matrix form. ⎡ Xi ,1 ⎤ ⎡ U1,1 U2,1 ⎢ ⎥ ⎢ ⎢ Xi ,2 ⎥ ⎢ U1,2 U2,2 ⎢ ⎥=⎢ ⋮ ⎢ ⋮ ⎥ ⎢ ⋮ ⎢ X ⎥ ⎢U ⎣ i , Nu ⎦ ⎣ 1, Nu U2, Nu

⋯ UNu ,1 ⎤⎡ Gi ,1 ⎤ ⎥⎢ ⎥ ⋯ UNu ,2 ⎥⎢ Gi ,2 ⎥ ⎥⎢ ⎥ ⋱ ⋮ ⎥⎢ ⋮ ⎥ ⎥ ⋯ UNu , Nu ⎦⎢⎣Gi , Nu ⎥⎦

(4)

Correlated inputs of the following form can be applied for the identification scheme. Uj , l = αj , lUl

⎡ Xi ,1 ⎤ ⎡ α1,1U1 α2,1U1 ⎢ ⎥ ⎢ ⎢ Xi ,2 ⎥ ⎢ α1,2U2 α2,2U2 ⎢ ⎥=⎢ ⋮ ⎢ ⋮ ⎥ ⎢ ⋮ ⎢ X ⎥ ⎢α U α U ⎣ i , Nu ⎦ ⎣ 1, Nu Nu 2, Nu Nu

1, 2, ..., Ny are the noise-free outputs; and yi are the corresponding measured outputs. The dimension of the MIMO system is Nu × Ny. The basic input−output relations for the system can be described using the following Laplace domain expressions.

∑ Gi ,j(s)Uj(s)

Yi(s) = Xi(s) + Ei(s)

⎡ Xi ,1 ⎤ ⎡U1 ⎥ ⎢ ⎢ ⎢ Xi ,2 ⎥ ⎢ 0 ⎥=⎢ ⎢ ⎢ ⋮ ⎥ ⎢⋮ ⎢X ⎥ ⎢ 0 ⎣ i , Nu ⎦ ⎣

(1)

i = 1, 2, ..., Ny

α Nu ,1U1 ⎤⎡ Gi ,1 ⎤ ⎥⎢ ⎥ ⋯ α Nu ,2U2 ⎥⎢ Gi ,2 ⎥ ⎥⎢ ⎥ ⋱ ⋮ ⎥⎢ ⋮ ⎥ ⎥ ⋯ α Nu , NuUNu ⎦⎢⎣Gi , Nu ⎥⎦ ⋯

(5)

The input matrix on the right side can be decomposed into a diagonal matrix containing the inputs and a coefficient matrix containing the constant multiplier terms.

Nu j=1

l = 1, 2, ..., Nu

So for the lth experiment, the input for any channel is a scaled form of a common input signal Ul. Using similar notations for all experiments the following equation is obtained.

Figure 2. Block diagram of a MIMO open-loop process with Nu inputs and Ny outputs.

X i (s ) =

j = 1, 2, ..., Nu

(2)

where Gi,j is the transfer function between the ith output and jth input. Here, upper case letters are used for variables in the Laplace (s) domain, i.e., U(s), X(s), and Y(s) denote the input, noise-free output, and noisy measured output, respectively. The corresponding lower case letters are used for variables in the time (t) domain. However, the indices (s) and (t) denoting the two domains will be frequently omitted for simplicity. Ei represents the noise term, and it is assumed that the measurement noise as well as the disturbances in an output signal that are not from the input excitation can be lumped into one additive term Ei. For a process with Nu inputs, Nu number of different experiments are required to perform in this approach. One or more or all of the inputs can be changed at each experiment. Irrespective of how the inputs are changed, the mathematical formulation will be presented based on the assumption that all the inputs are changed during each experiment. So for the lth experiment, the process is excited through its different input channels using the input signals uj,l with j = 1, 2, ..., Nu and the output measurements are recorded. Nu similar experiments are to be carried out using different sets of inputs for each experiment. The input−output relations for the lth experiment can be described by

0 ⎤⎡ α1,1 α2,1 ⎥⎢ U2 ⋯ 0 ⎥⎢ α1,2 α2,2 ⎥⎢ ⋮ ⋱ ⋮ ⎥⎢ ⋮ ⋮ 0 ⋯ UNu ⎥⎦⎢⎣ α1, Nu α2, Nu 0 ⋯

⋯ α Nu ,1 ⎤⎡ Gi ,1 ⎤ ⎥ ⎥⎢ ⋯ α Nu ,2 ⎥⎢ Gi ,2 ⎥ ⎥ ⎥⎢ ⋱ ⋮ ⎥⎢ ⋮ ⎥ ⋯ α Nu , Nu ⎥⎦⎢Gi , N ⎥ ⎣ u⎦ (6)

Note that in the coefficient matrix, if the off-diagonal terms are set to zeros, the experiments will be equivalent to sequential tests where inputs are changed one at a time. Next, the above equation is premultiplied by the inverse of the input data matrix. The multiplication results in the expression of the transfer functions separated from the data matrix. ⎡ Xi ,1 ⎤ ⎢ ⎥ ⎢ U1 ⎥ ⎢ X ⎥ ⎡ α1,1 α2,1 ⎢ i ,2 ⎥ ⎢ α α2,2 ⎢ U2 ⎥ = ⎢ 1,2 ⎢ ⎢ ⎥ ⋮ ⋮ ⎢ ⋮ ⎥ ⎢ α α ⎢ ⎢ X ⎥ ⎣ 1, Nu 2, Nu ⎢ i , Nu ⎥ ⎢⎣ UN ⎥⎦ u

⋯ α Nu ,1 ⎤⎡ Gi ,1 ⎤ ⎥ ⎥⎢ ⋯ α Nu ,2 ⎥⎢ Gi ,2 ⎥ ⎥ ⎥⎢ ⋱ ⋮ ⎥⎢ ⋮ ⎥ ⋯ α Nu , Nu ⎥⎦⎢Gi , N ⎥ ⎣ u⎦ (7)

In another step of premultiplication the equation is multiplied by the inverse of the matrix containing the input coefficients. C


Article

Industrial & Engineering Chemistry Research ⎡ α1,1 α2,1 ⎢ ⎢ α1,2 α2,2 ⎢ ⋮ ⎢ ⋮ ⎢α α ⎣ 2, Nu 2, Nu

⋯ α Nu ,1 ⎤ ⎥ ⋯ α Nu ,2 ⎥ ⎥ ⋱ ⋮ ⎥ ⋯ α Nu , Nu ⎥ ⎦

−1

⎡ α1,1 α1,2 ̅ ̅ ⎢ α α ⎢ ̅2,1 ̅2,2 =⎢ ⋮ ⎢ ⋮ ⎢α ⎣ ̅Nu ,1 α̅Nu ,2

⋯ α1,̅ Nu ⎤ ⎥ ⋯ α̅2, Nu ⎥ ⎥ ⋱ ⋮ ⎥ ⋯ α̅Nu , Nu ⎥⎦

⎡ α1,1 α1,2 ̅ ̅ ⎢ α α ⎢ ̅2,1 ̅2,2 ⎢ ⋮ ⎢ ⋮ ⎢α ⎣ ̅Nu ,1 α̅Nu ,2

(8)

Comment: It is required that the coefficients matrix should be invertible. This implies that the identification experiments should be unique, i.e., the inputs should be changed differently at each experiment. This step of multiplication expresses the transfer functions relating the ith output to all the inputs in terms of the input− output data.

⎡ G1,1 G2,1 ⎢ ⎢ G1,2 G2,2 ⎢ ⎢ ⋮ ⋮ ⎢ ⎢⎣G1, Nu G2, Nu

Nu

∑ αj̅ ,l l=1

X i , l (s ) (11)

⋯ ⋯ ⋱ ⋯

⎡X X 2,1 ⎢ 1,1 U1 ⎢ U α1,̅ Nu ⎤⎢ 1 ⎥⎢ X X 2,2 α̅2, Nu ⎥⎢ 1,2 U2 ⎥ U2 ⋮ ⎥⎢⎢ ⋮ ⋮ α̅Nu , Nu ⎥⎦⎢ ⎢ X1, Nu X 2, Nu ⎢ ⎣ UNu UNu

Nu

hl s

Y (s ) =

(12)

where hl is the size of the step. So for the step inputs 1 Gi , j(s) = s

∑ l=1

αj̅ , l hl

Nu

E (s ) =

X i , l (s )

1 = s

Nu

∑ l=1

αj̅ , l

(13)

hl

[Yi , l(s) − Ei , l(s)]

bms m + bm − 1s m − 1 + ⋯+ b0 n

s + an − 1s

n−1

+ ⋯+ a0

X Ny ,1 ⎤ ⎥ U1 ⎥ ⎥ X Ny ,2 ⎥ ⋯ ⎥ U2 ⎥ ⋱ ⋮ ⎥ ⎥ X Ny , Nu ⎥ ⋯ ⎥ UNu ⎦ ⋯

(10)

αj̅ , l hl αj̅ , l hl

Yi , l(s)

(16)

Ei , l(s)

(17)

Note that it is not required to evaluate Y(s); instead the parameter estimation equation will be formulated in terms of y(t) that can be evaluated from the time domain measurements of the ith output signal. Using the above notations a generalized parameter estimation equation can be written as follows. Y (s ) =

bms m + bm − 1s m − 1 + ⋯+ b0 n

s + an − 1s

n−1

+ ⋯+ a0

e−δs

1 + E (s ) s

(18)

In the equation error format with V(s) being the equation error term, the above equation can be rearranged to give

(14)

[s n + an − 1s n − 1 + ⋯+ a0]Y (s)

The above equation formulates the basis for the parameter estimation of the transfer function Gi,j(s) which can be expressed as Gi , j(s) =

∑ l=1

Note that αj̅ , l and hl being known constants, the right-hand side term of the above equation is the linear combination of the ith output signals for different experiments. As seen from the lefthand side of the expression, the linear combination should generate the unit step response for the corresponding transfer function. However, Xi,l are not available, rather we have the output measurements as Yi,l. If we assume that the noise terms in Yi,l can be combined to a single term Ei,l, we have Gi , j(s)

∑ l=1

Nu

(9)

Here, ak, k = 0, ···, n are the denominator coefficients, bp, p = 0, ···, m are the numerator coefficients, and δ is the time delay. Without loss of generality, the coefficient an can be normalized to be 1. No reference to the indices (i,j) is made for parameters a, b, and δ in the above expression for simplicity in the presentation. Also the following notations are defined for the sake of simplicity in the presentation of the parameter estimation equation.

For step input signals we get Ul =

⋱ ⋯

for all of the outputs giving the following expression.

⋯ G Ny ,1 ⎤ ⎡ α1,1 α1,2 ̅ ⎥ ⎢ ̅ ⎥ ⋯ G Ny ,2 α̅2,2 ⎢α ⎥ = ⎢ ̅2,1 ⋮ ⋱ ⋮ ⎥ ⎢ ⋮ ⎥ ⎢ α α ⋯ G Ny , Nu ⎥⎦ ⎣ ̅Nu ,1 ̅Nu ,2

Ul(s)

⋯

Following the same procedure, similar equations can be derived

From the above equation the individual transfer functions can be expressed as the following. Gi , j(s) =

⋯

⎡ Xi ,1 ⎤ ⎢ ⎥ U α1,̅ Nu ⎤⎢ 1 ⎥ ⎡ Gi ,1 ⎤ ⎥ ⎥⎢ Xi ,2 ⎥ ⎢ ⎥ ⎢G ⎥ α̅2, Nu ⎥⎢ i ,2 ⎥ ⎥⎢ U2 ⎥ = ⎢ ⎥ ⎢ ⋮ ⎥ ⋮ ⎥⎢ ⎢ ⋮ ⎥ ⎢ α̅Nu , Nu ⎥⎦⎢ G ⎥ Xi , Nu ⎥ ⎣ i , Nu ⎦ ⎢ ⎥ ⎢⎣ UN ⎥⎦ u

= [bms m + bm − 1s m − 1 + ⋯+ b0]e−δs

1 + V (s ) s

(19)

Equation 19 is the general form for parameter estimation of SISO transfer function models in the equation error format. A number of approaches is available to estimate the model

e−δs (15) D


Article

Industrial & Engineering Chemistry Research an − 1 ⎡ ⎤ ⎢ ⎥ ⋮ ⎢ ⎥ ⎢ ⎥ a0 ⎢ ⎥ b0 ⎢ ⎥ ⎢ min(m ,1) ⎥ ⎢ ( −δ)1 − q ⎥ ⎢ ∑ bq ⎥ (1 − q)! ⎥ θ = ⎢ q=0 ∈ (2n + 1) × 1 ⎢ ⎥ ⎢ min(m ,2) ( −δ)2 − q ⎥ ⎢ ∑ bq ⎥ (2 − q)! ⎥ ⎢ q=0 ⎢ ⎥ ⋮ ⎢ ⎥ ⎢ m n−q ⎥ ⎢ ∑ bq ( −δ) ⎥ ⎢⎣ q = 0 (n − q)! ⎥⎦

parameters. We will follow the integral equation approach37,38 as it allows to simultaneously estimate the time delay along with the numerator and denominator coefficients. Following the integral equation approach, dividing both sides of the equation by sn we get n

Y (s ) +

∑ an− k k=1

m

Y (s ) s

=

k

∑ bm− p p=0

1 s

n−m+p+1

e−δs + Vn(s) (20)

where Vn(s) = V(s)/sn. Note that division by sn in Laplace domain is equivalent to nth order integration of the time domain equation. Equation 20 can be expressed in the time domain as n

y(t ) +

m

∑ an− ky[k](t ) = ∑ bm− p k=1

p=0

Equation 26 can be written for t = td+1, td+2, ···, tN and combined to give

(t − δ)n − m + p + vn(t ) (n − m + p )!

Γ = Φθ + Vn

(21)

⎡ γ (t ) ⎤ ⎢ d+1 ⎥ ⎢ γ(td + 2)⎥ Γ(t ) = ⎢ ⎥, ⎢ ⋮ ⎥ ⎢ ⎥ ⎣ γ(tN ) ⎦

t

∫0

y[k](t ) =

∫0 ∫0

y(τ ) dτ t

τ1

⋯

∫0

(22) τk − 1

y(τk) dτk ⋯ dτ1

for

k≥2 (23)

Nu

∑

αj̅ , l

l=1

hl

yi , l (t )

(24)

Following the mathematical manipulation in Ahmed et al.37 we

θ LS = (ΦT Φ)−1ΦT Γ

can write n

y(t ) = − ∑ an − k y[k](t ) + b0 k=1

tn + n!

n

∑ p=1

⎛ min(m , p) p−q ⎞ ⎜ ∑ b ( − δ ) ⎟ + v (t ) q n ⎜ (p − q)! ⎟⎠ ⎝ q=0

(25)

(26)

where γ(t ) = y(t )

(27)

⎤ ⎡ tn t n−1 ϕ(t ) = ⎢−y[1](t ) ⋯ −y[n](t ) ⋯ t 1⎥ n! (n − 1)! ⎦ ⎣ ∈ 1 × (2n + 1)

(31)

(32)

The error terms in the estimation equation evolve due to the presence of noise in the output measurements. Typically, the measurement noise is zero mean white noise or filtered white noise. In the presence of colored noise, the least-squares solution may not be unbiased. Even if the measurement noise is assumed to be white with zero-mean, the integration operation performed on the output signal results in a colored error term. So, the LS solution may not be unbiased even for a white measurement noise, and we need a bias elimination scheme. To get an unbiased estimate of the parameters, different techniques can be used. We use the instrumental variable (IV) method which is commonly used in continuous-time identification.37,39 To generate the instruments, parameters obtained using the LS solution are used to get the predicted values of the output. The instrument matrix is then derived by replacing the terms related to the output, y(t), in the regressor by their predicted values, y ̂(t ), which is the predicted unit step response of the (i,j)th transfer function. The response is predicted using the estimated parameters and delay obtained from the least-squares solution. For the nth order process, the instrument vector is defined as

t n−p (n − p)!

Or equivalently γ(t ) = ϕ(t )θ + vn(t )

⎡ ϕ(t )⎤ ⎢ d+1 ⎥ ⎢ ϕ(td + 2)⎥ Φ(t ) = ⎢ ⎥ ⎢ ⋮ ⎥ ⎢ ⎥ ⎣ ϕ(tN ) ⎦

Here, d is the time delay in terms of number of sampling intervals (Δt), i.e., d = δ/Δt and N is the total number of samples available. When the time delay is not an integer multiple of the sampling interval, d is chosen as the nearest integer in the positive direction. The least-squares (LS) solution for the estimation eq 30 is given by

Also from eq 16 y(t ) =

(30)

with

where y[k](t) denotes the kth order integration of y(t), i.e., y[1](t ) =

(29)

(28) E


Article

Industrial & Engineering Chemistry Research ⎤ ⎡ tn t n−1 ψ (t ) = ⎢−y[1] ⋯ t 1⎥ ̂ (t ) ⋯ −y[̂ n] (t ) n! (n − 1)! ⎦ ⎣

manipulated by changing its fuel feed. Both of the inputs are measured in terms of the % openings of the corresponding valves while the compositions are measured as the mol % of benzene. Figure 3 shows the interface for simulation of the column in Loop Pro.

(33)

An alternative formulation of the instrument matrix can be found in ref 38. The instrumental variable estimate of the parameters is given by θ IV = (ΨT Φ)−1ΨT Γ

(34)

The IV procedure can be repeated. From the solution of eq 34, an−1, ···, a0 can be obtained directly from the first n elements of θn as an−k = θn(k), k = 1, ···, n. b0 is also obtained directly as b0 = θn+1. For a first order process with n = 1 and m = 0, the delay δ can be obtained as δ = −θ3/θ2

(35)

For processes with n = 2 and m = 1 δ=

−θ4 ±

θ4 2 − 2θ3θ5 θ3

b1 = θ4 + θ3δ

(36) (37)

Figure 3. Simulation interface of the distillation column in Loop Pro.

For higher order processes to extract δ the following polynomial is required to solve. m+1

∑ i=0

δ m+1−i θn + 1 + i = 0 ( m + 1 − i )!

(38)

The estimated δ can then be used to get the estimate of b1, ···, bm using ⎡ ⋯ 0 ⎤⎡ 1 ⎡ b1 ⎤ ⎢ δ ⎥ θn + 1 ⎤ ⋮ ⋱ ⋮ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⋮ ⎢ ⋮ ⎥ ⎢⋮⎥ = ⎢ m ⎥ m−1 δ ⎢ ⎥ ⎢δ ⋯ 1 ⎥⎢⎣ θn + m + 1⎥⎦ ⎣bm ⎦ ⎢ ⎥⎦ ! − ! m ( m 1) ⎣

(39)

Thus, all of the parameters of the transfer function Gi,j(s) are estimated. The same procedure is followed for all of the transfer functions for i = 1, 2, ..., Ny and j = 1, 2, ..., Nu, resulting the MIMO transfer function matrix.

3. APPLICATION TO A SIMULATED DISTILLATION COLUMN The proposed methodology was used to identify a MIMO model of a simulated distillation column. In this section, the identification results and the closed loop response of the process under a dynamic matrix controller (DMC) using the identified model are presented. The distillation column is available as a simulation module in the Loop Pro software [Loop Pro is a training simulator, developed and marketed by the Control Station Inc., U.S.A. (www.controlstation.com)]. A binary solution of benzene and toluene in considered in this system. The objective of the operation is to send a high percentage of benzene (and thus low percentage of toluene) out the tops stream and a low percentage of benzene (and thus high percentage of toluene) out the bottoms stream. The column dynamic model employs tray-by-tray mass and energy calculations similar to that proposed in ref 40. The outputs considered for this study are the composition of the top product, yD, and that of the bottom product, xB. The inputs are the reflux rate, R, and the reboiler duty, QB. The reflux rate is manipulated by the reflux valve while the reboiler duty is

A MIMO model of the distillation column, consisting of four transfer functions, was estimated using the proposed methodology by simultaneously changing both of the inputs. As there are two inputs, a set of two experiments was required to be performed. The step sizes for the two experiments were [R QB] = [2 1] and [2 −1]. So in terms of the notations used in this paper U1 =

1 s

U2 =

1 s

⎡ α1,1 α2,1 ⎤ ⎡ 2 1 ⎤ ⎢ ⎥=⎢ ⎥ ⎣ α1,2 α2,2 ⎦ ⎣ 2 −1⎦ (40)

eq 41 shows the transfer functions estimated using the proposed methodology. ⎡ 0.8187e−22.6s −0.6053e−16s ⎤ ⎢ ⎥ ⎡ yD ⎤ ⎢ 71s + 1 61.3s + 1 ⎥⎡ R ⎤ = ⎢ ⎥ ⎢ ⎥⎢ ⎥ −27.6s ⎣ xB ⎦ −0.6887e−17.9s ⎥⎣Q B ⎦ ⎢ 0.6e ⎢⎣ 85.2s + 1 95.8s + 1 ⎥⎦ F

(41)


Article


Figure 4. Closed loop response of the distillation column under a DMC controller using the two sets of models. Top-left: response in yD; top-right: response in xB; bottom-left: input move in R; bottom-right: input move in QB.

(eq 41) is used in the DMC controller. Case 2 refers to the case when the model estimated by sequentially changing one input at a time (eq 42 is used). The responses show that a better tracking performance is obtained in case 1 than in case 2. Relatively more aggressive control action is observed for case 2 resulting in larger overshoot of both of the outputs compared to case 1.

To demonstrate the efficacy of the proposed algorithm, transfer functions for the input−output pairs were also estimated by changing one variable at a time. Equation 42 shows the estimated transfer functions. In this case, first R was changed by 2% and the transfer functions between yD and R and between xB and R were estimated. Next, R was kept constant while QB was changed. As for the previous experiment, two experiments were conducted by changing QB by +1% and −1%, in this case QB was also changed by +1% and −1% to make the input excitation equivalent for the two cases. The transfer functions between QB and the two outputs are the averages of the transfer functions obtained using the positive and the negative changes. ⎡ 0.8548e−25s −0.83e−19s ⎤ ⎢ ⎥ ⎡ yD ⎤ ⎢ 74.5s + 1 70.4s + 1 ⎥⎡ R ⎤ ⎢ ⎥=⎢ ⎥⎢ ⎥ −28.4s ⎣ xB ⎦ −0.39e−11s ⎥⎣Q B ⎦ ⎢ 0.5757e ⎣ 88s + 1 67.5s + 1 ⎦

4. SIMULATION RESULTS To demonstrate the performance of the proposed methodology in terms of the properties of the estimated parameters and the robustness of the algorithm, a set of simulation results is presented. Numerical calculations are carried out using Matlab and Simulink. In Matlab the lsim command is used with the system defined as continuous-time transfer functions. In Simulink processes are defined as continuous-time transfer function with the solver chosen as ODE45. Noise free data were first generated, and white noise sequences were added to the noise free signal to simulate noisy scenarios. For the following results, noise to signal ratios were 10%. Monte Carlo simulations (MCS) were performed to obtain the properties of the estimated parameters. This is done in Matlab by changing the seed of the added noise. The step responses were continued until a final steady state is reached. For each step response 500 data points were collected. The sampling intervals were between 0.1 and 0.2 to meet the requirements of data length. 4.1. Identification of a 3 × 3 System. A 3 × 3 system is used to demonstrate the properties of the estimated parameters. Equation 43 shows the transfer function matrix of the process.

(42)

A considerable difference in the values of the model parameters estimated using the two approaches is noticed; especially for the two models related to the reboiler duty. To compare the performance of the models identified using the proposed methodology and that obtained by changing one variable at a time, a DMC controller is implemented to the process using the identified models. Figure 4 shows the closed loop responses in yD and xB under a DMC controller. The set point was changed simultaneously from 94.5% to 95.5% for yD and that from 2.6% to 2.1% for xB. Case 1 refers to the case when the model estimated using the proposed methodology

⎡ 1.986e−0.71s − 5.24e−6s − 5.984e−2.24s ⎢ 26.7 s + 1 34.9 s + 1 14.29s + 1 ⎢ ⎢ −0.68s −3.59s (− 1s + 0.33)e − 2.38e−0.42s ⎢ − 0.0204e ⎢ 50.98s 2 + 10s + 1 40.7s 2 + 2.55s + 1 204.5s 2 + 22.9s + ⎢ ⎢ − 0.374e−7.75s − 11.3e−3.79s 9.811e−1.59s ⎢⎣ 2 22.22s + 1 11.36s + 1 472.6s + 43.48s + 1

G

⎤ ⎥ ⎥ ⎥ ⎥ 1⎥ ⎥ ⎥ ⎥⎦

(43)


Article


Figure 5. A typical data set from the identification experiment for the 3 × 3 system.

Figure 6. Nyquist plots of the different transfer functions. + indicates the plots of the true process and the lines represents the 100 estimated models.

rows represent the three outputs. Equation 44 shows the identified models. The results are from 100 MCS. The mean values of the parameter estimates are presented along with the standard deviations of the parameters in the parentheses.

Figure 5 shows a typical data set required for identification of the nine transfer functions for the 3 × 3 system. The individual columns in the figure represents individual experiments. As there were 3 inputs, 3 sets of experiments were performed. The

⎡ 1.986(± 0.001)e−0.71(±0.028)s − 5.25(± 0.07)e−6(±0.6)s ⎢ 26.7(± 0.06)s + 1 34.99(± 1.76)s + 1 ⎢ ⎢ −3.59(±0.27)s ⎢ − 0.0204(± 0.000)e (− 0.98(± 0.12)s + 0.33(± 0.001))e−0.82(±0.5)s ⎢ 2 40.5(± 1.2)s 2 + 2.53(± 0.1)s + 1 ⎢ 50.8(± 3.16)s + 9.97(± 0.2)s + 1 ⎢ −7.77(±0.17)s − 11.3(± 0.23)e−4.05(±1.32)s ⎢ − 0.374(± 0.001)e ⎢⎣ 22.19(± 0.33)s + 1 460(±73)s 2 + 43.4(± 1.6)s + 1

4.2. Identification of a 4 × 4 System. A 4 × 4 process is used to demonstrate the frequency domain properties of the identified models. Figure 6 shows the Nyquist plots of the true transfer functions and that of the 100 identified transfer functions for each input−output pair from 100 Monte Carlo simulations. Transfer functions used to generate data were

⎤ ⎥ ⎥ ⎥ −0.42(±0.14)s ⎥ − 2.38(± 0.002)e ⎥ 204.46(± 3.81)s 2 + 22.88(± 0.1)s + 1 ⎥ ⎥ 9.81(± 0.03)e−1.62(±0.2)s ⎥ ⎥⎦ 11.32(± 0.31)s + 1 − 5.988(± 0.03)e−2.25(±0.29)s 14.32(± 0.46)s + 1

(44)

either first order or overdamped second order. All of the estimated transfer functions were first order. The overlapping of the frequency response of the estimated models with one another and with that of the true process shows the accuracy and consistency in the parameter estimates. The mismatches H


Article


Figure 7. Estimated parameters of the four transfer functions of the Wood and Berry distillation column; clock-wise from top left yD − R, yD − QB, xB − R, xB − QB.

for the two models Y3−U1 and Y3−U2 resulted due to the first order approximation of the second order processes. 4.2.1. Effect of Noise. To demonstrate the robustness of the proposed method, identification results for different levels of noise are presented in this section. An often cited model of a distillation column, namely, the Wood and Berry column,41 having the following transfer function is used for this study. ⎡ 12.8e−1s −18.9e−3s ⎤ ⎢ ⎥ ⎡ yD ⎤ ⎢ 16.7s + 1 21s + 1 ⎥⎡ R ⎤ ⎢ ⎥=⎢ ⎥⎢ ⎥ − 7s ⎣ xB ⎦ −19.4e−3s ⎥⎣Q B ⎦ ⎢ 6.6e ⎣ 10.9s + 1 14.4s + 1 ⎦

5. CONCLUDING REMARKS MIMO models are often identified by changing the inputs sequentially. While this simplifies the parameter estimation method, the identified models result in poor performance when used for control. This article proposes a new methodology to identify MIMO models that allows all the inputs to be changed simultaneously. At the same time, parameters of each transfer function are estimated separately by solving independent equations. Control performance of the identified models has been demonstrated using the example of a simulated distillation column under a model predictive controller. The accuracy and consistency in parameter estimation is demonstrated using a set of simulation results. Future efforts will focus on extending the developed methodology for other types of input signals.

(45)

■

Figure 7 shows the identification results for different levels of noise added to the output signals. The individuals plots are for each of the four transfer functions. The presented parameters are the mean values of 100 Monte Carlo simulation estimates. The vertical error bars show the standard deviations of the 100 estimates for the corresponding parameters. The dotted lines indicate their true values. The simulation results show that for noise level as high as NSR = 25% satisfactory results are obtained. The quality of parameter estimation may be the result of data quality or data processing. Higher sampling time with high NSR may make the estimates poor. The IV algorithm was repeated a fixed 5 times for each case as there was no significant improvement beyond that. In this article the main contribution is in the conversion of the MIMO identification problem into a SISO parameter estimation problem. Consequently identification results of MIMO systems with different dimensions are presented for low order transfer functions. Study of the effect of sampling interval has been done in the literature for similar methods. Although the examples only cover square systems meaning that the number of inputs equals the number of outputs, the methodology is applicable for nonsquare systems as well.

AUTHOR INFORMATION

Corresponding Author

*(S. Ahmed) Phone: +1 709 864 7652. Fax: +1 709 864 7858. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the financial support from Research and Development Corporation (RDC) of Newfoundland and Labrador and Natural Sciences and Engineering Research Council (NSERC) of Canada.

■

REFERENCES

(1) Waller, J.; Waller, K. Defining directionality: Use of directionality measures with respect to scaling. Ind. Eng. Chem. Res. 1995, 34, 1244− 1252. (2) Zhu, Y.; Stec, P. Simple control-relevant identification test methods for a class of ill-conditioned processes. J. Process Control 2006, 16, 1113−1120. I


Article


continuous-time models from sampled data; Garnier, H., Wang, L., Eds.; Springer: London, 2008; pp 339−362. (27) Kim, Y. C.; Jin, L. Robust identification of continuous-time loworder models using moments of a single rectangular pulse response. J. Process Control 2013, 23, 682−695. (28) Doma, M. J.; Taylor, P. A.; Vermeer, P. J. Closed loop identification of MPC models for MIMO processes using genetic algorithm and dithering one variable at a time: Application to an industrial distillation tower. Comput. Chem. Eng. 1996, 20, S1035− S1040. (29) Choi, J.; Lee, J.; Jung, J.; Lee, M.; Han, C. Sequential loop closing identification of multivariable process models. Comput. Chem. Eng. 2000, 24, 809−814. (30) Choi, J.; Lee, J.; Edgar, T. Use of the sequential loop closing method for iterative identification of ill-conditioned processes. Ind. Eng. Chem. Res. 2000, 39, 2404−2409. (31) Li, S. Y.; Cai, W. J.; Mei, H.; Xiong, Q. Robust decentralized parameter identification for two-input two-output processes from closed-loop step responses. Control Engineering Practice 2005, 13, 519−531. (32) Dayal, B.; McGregor, J. Multi-input process identification. J. Process Control 1997, 7, 269−282. (33) Conner, J.; Seborg, D. An evaluation of MIMO input designs for process identification. Ind. Eng. Chem. Res. 2004, 43, 3847−3854. (34) Jin, Q.; Wang, Z.; Yang, R.; Wang, J. An effective direct closed loop identification method for linear multivariable systems with colored noise. J. Process Control 2014, 24, 485−492. (35) Sun, H.; Betti, R. Simultaneous identification of structural parameters and dynamic input with incomplete output-only measurements. Structural Control and Health Monitoring 2014, 21, 868−889. (36) Sun, H.; Lus, H.; Betti, R. Identification of structural models using a modified Artificial Bee Colony algorithm. Comput. Struct. 2013, 116, 59−74. (37) Ahmed, S.; Huang, B.; Shah, S. L. Identification from step response with transient initial conditions. J. Process Control 2008, 18, 121−130. (38) Wang, Q. G.; Zhang, Y. Robust identification of continuous systems with dead-time from step responses. Automatica 2001, 37, 377−390. (39) Ahmed, S.; Huang, B.; Shah, S. L. Novel identification method from step response. Control Engineering Practice 2007, 15, 545−556. (40) McCune, L. C.; Gallier, P. W. Digital Simulation: A Tool for Analysis and Design of Distillation Column. ISA Trans. 1973, 12, 193. (41) Wood, R. K.; Berry, M. W. Terminal composition control of a binary distillation column. Chem. Eng. Sci. 1973, 28, 1707−1717.

(3) Skogestad, S.; Morari, M.; Doyle, J. Robust Control of IllConditioned Plants: High-Purity Distillation. IEEE Trans. Autom. Control 1988, 33, 1092−1105. (4) Skogestad, S.; Lundstrom, P.; Jacobson, E. W. Selecting the best distillation column configuration. AIChE J. 1990, 36, 753−764. (5) Luyben, W. Derivation of transfer functions for highly nonlinear distillation columns. Ind. Eng. Chem. Res. 1987, 26, 2490−2495. (6) Martin, P.; Odloak, D.; Kassab, F. Robust model predictive control of a pilot plant distillation column. Control Engineering Practice 2013, 21, 231−241. (7) Koung, C. W.; MacGregor, J. F. Design of Identification Experiments for Robust Control. A Geometric Approach for Bivariate Processes. Ind. Eng. Chem. Res. 1993, 32, 1658−1666. (8) Rao, G. P.; Garnier, H. Numerical illustrations of the relevance of direct continuous-time model identification. Proceedings of the 15th IFAC Triennial Congress, Barcelona, Spain, July 2002; IFAC: 2002. (9) Sung, S.; Lee, I. Prediction error identification method for continuous-time processes with time delay. Ind. Eng. Chem. Res. 2001, 40, 5743−5751. (10) Bastogne, T.; Garnier, H.; Sibille, P. A PMF-based subspace method for continuous-time model identification. Application to a multivariable winding process. Int. J. Control 2001, 74, 118−132. (11) Johansson, R.; Verhaegen, M.; Chou, C. T. Stochastic theory of continuous-time state-space identification. IEEE Transaction on Signal Processing 1999, 47, 41−51. (12) Li, W.; Raghavan, H.; Shah, S. L. Subspace identification of continuous-time models for process fault detection and isolation. J. Process Control 2003, 13, 407−421. (13) Ohsumi, A.; Kameyama, K.; Yamaguchi, K. Subspace identification for continuous-time stochastic systems via distribution based approach. Automatica 2002, 38, 63−79. (14) Aziz, M.; Mohd-Mokhtar, R.; Wang, L. Identification step response estimates utilizing continuous time subspace approach. J. Process Control 2013, 23, 254−270. (15) Jeng, J.; Lin, Y. Closed-loop identification of dynamic models for multivariate systems with applications to monitoring and redesign of controllers. Ind. Eng. Chem. Res. 2011, 50, 1460−1472. (16) Melo, D. L.; Friedly, J. C. On-line, Closed-loop identification of multivariable systems. Ind. Eng. Chem. Res. 1992, 31, 274−281. (17) Mei, H.; Li, S.; Cai, W.-J.; Xiong, Q. Decentralized closed-loop parameter identification for multivariable processes from step responses. Mathematics and Computer in Simulation 2005, 68, 171− 192. (18) Mei, H.; Li, S. Decentralized identification of multivariable integrating processes with time delay from closed-loop step tests. ISA Trans. 2007, 46, 189−198. (19) Wang, Q. G.; Zhang, Y. A novel FFT-based robust multivariate process identification method. Ind. Eng. Chem. Res. 2001, 40, 2485− 2494. (20) Zhu, Y. Parametric versus nonparametric models in MPC process identification. Hydrocarbon Processing 2000, 79, 65−72. (21) Jin, Q.; Cheng, Z.; Dou, J.; Cao, L.; Wang, K. A novel closed loop identification method and its application of multivariable system. J. Process Control 2012, 22, 132−144. (22) Rajapandiyan, C.; Chidambaram, M. Closed-loop identification of multivariable systems by optimization method. Ind. Eng. Chem. Res. 2012, 51, 1324−1336. (23) Rajapandiyan, C.; Chidambaram, M. Closed-loop identification of second order plus time delay model of multivariable systems by optimization method. Ind. Eng. Chem. Res. 2012, 51, 9620−9633. (24) Jeng, J. Closed-loop identification of multivariable systems using B-Spline series expansion for step response. Ind. Eng. Chem. Res. 2012, 51, 2376−2387. (25) Wang, Y.; Cai, W.; Ge, M. Decentralized relay-based multivariable process identification in the frequency domain. IEEE Trans. Autom. Control 2003, 48, 873−877. (26) Yang, Z. Iterative methods for identification of multiple-inputt continuous-time systems with unknown time delays. In Identification of J


Identification of MIMO Continuous-Time Models Using Simultaneous

Recommend Documents