On closed loop transient response system identification

Some methods for transient closed loop step response system identification presented in the literature are reviewed. Interestingly some errors in a method published in the early 80’s where propagated into a recently published method. These methods are reviewed and some improved methods are suggested and presented. The methods are compared against each other on some closed loop system examples, e.g. a well pipeline-riser severe-slugging flow regime example, using Monte Carlo simulations for comparison of the methods.


Introduction
Closed loop subspace System IDentification (SID) has been given a large amount of effort the last two decades resulting in e.g. the DSR e alg., Di Ruscio (2009). A problem with closed loop subspace SID is the possibility of correlations between the controller output and the process noise, if this is present then the results may be biased. However, this is solved in DSR e, see Di Ruscio (2009). Closed loop SID is a relatively matured field and is documented in books as e.g. Ljung (1998), Söderström and Stoica (1989) etc.
The open loop plant may be unstable or very poorly damped and has to be operated in closed loop, or it might be for economic or safety reasons (Söderström and Stoica (1989)).
In this paper our main focus is the SID method used in so-called on-line closed loop step response tuning methods: Yuwana and Seborg (1982), Jahanshahi and Skogestad (2015).
We can classify such methods as indirect, as we need to know the reference and output time series, as well as the controller (See Fig.1). We may further call such methods nonparametric, type transient analysis, as described in Söderström and Stoica (1989), which are examples of continuous identification methods.
The only notable strength of these methods is that they are easy to use, pragmatic and easy to understand.
However, the weaknesses are rather many. Firstly, nonparametric transient step response methods are sensitive to noise and do not give accurate results (Söderström and Stoica (1989) Ch.3). Secondly, one requires a high number of measurable quantities sampled at high frequency, e.g. a settled output time series response. Third, some of these curve points are rather strenuous to identify, e.g. the first valley point. Producing larger oscillations by increasing the reference change or the controller proportional gain could help, however this may be rather unappetizing in process industries.
One contribution of this paper is to show that the identification method given in Jahanshahi and Skogestad (2015) may be improved. It should be a well-known fact that if we are given an experiment signal using a h c (s) h p (s)

Method 1
Method 2 r u y − + Figure 1: Block diagram of a closed loop system which is illustrating two SID methods. Method 1: Direct method using y and u. Method 2: Indirect method using r, y and h c . Note that we also have a third approach called joint input-output method, where r, u and y are needed.
pure step signal only, then this signal is only persistent exciting of order one which is different from the second order model structure, assumed in Jahanshahi and Skogestad (2015). Furthermore, we may predict that by introducing process noise, this method will present itself inferior to other SID methods. Notice however that this particular SID method, i.e. Jahanshahi and Skogestad (2015), is aimed for a pipeline-riser system, where the severe-slugging flow regime is modeled using a second order marginally stable/unstable open loop model. Severe-slugging flow, i.e. oscillations in the outlet flow is regarded as a problem in the offshore industry and needs to be addressed for safety and economic reasons (Courbot (1996)).
The contributions of this paper can be itemized as follows • An error is discovered in Jahanshahi and Skogestad (2015) having it roots in Yuwana and Seborg (1982) for estimating the time instant of the first peak, t p , giving implications for the φ estimate in Jahanshahi and Skogestad (2015). The modified algorithm is denoted JSDR.
• An alternative formulation to the algorithm presented in Jahanshahi and Skogestad (2015), is proposed. The proposed algorithm is denoted DR.
• Documentation is given for an additional error discovered in Yuwana and Seborg (1982) for the time constant and the time delay expressions of a first order plus time delay model.
• Counter-example on the closed loop transient response SID algorithm presented in Jahanshahi and Skogestad (2015). Monte Carlo simulation experiments shows that the two transient response identification algorithms, JSDR and DR, considered in this paper, are inferior to the following subspace algorithms; DSR/DSR ry (Di Ruscio (1996)) and DSR e (Di Ruscio (2009)) and the Prediction Error Method (PEM).
• An m-file script developed for estimating the 6 necessary parameters from the transient closed loop step response e.g. t p etc, is presented in App. A.
The rest of the paper is organized as follows. In Sec. 2 we formulate the improved JSDR method, as well as the proposed alternative, the DR method. In Sec. 3 we give numerical examples. In Secs. 4 and 5 discussion and concluding remarks are given, respectively.

On transient response SID algorithms 2.1. Improved JSDR method
Assume that the underlying system can be described by the following transfer function This particular model was used to describe the essential dynamics of a pipeline-riser system in Jahanshahi and Skogestad (2015). By stabilizing with a P-controller, the closed loop transfer function becomes where we introduce a reference step, We restrict ourselves to 0 < ξ < 1 and τ > 0. The problem renders itself in identifying the closed loop model parameters (ξ, τ, K, τ z ), thereafter, by back calculation, we obtain the open loop model parameters (a 1 , a 0 , b 1 , b 0 ). These models are identified based on the data vectors, Figure 2: Closed loop step response. The x-axis t is time.
The estimated parameters (t p , y p , r, t u , y u , y ∞ ) are used to identify (ξ, τ, K, τ z ) in Eq. (2). A MATLAB function for identifying these six parameters is found in App. A. x-axis t is time. The estimated parameters (t p , y p , r, t u , y u , y ∞ ) are used to identify (ξ, τ, K, τ z ) in Eq.
(2). A MATLAB function for identifying these six parameters is found in App. A. generated from a closed loop step response as illustrated in Figs. 2 and 3. The time transient closed loop step response is where and In Jahanshahi and Skogestad (2015) the time response, Eq. (A.2), is wrong, even though referring to the correct Eq. (5) or Eq. (A-1) in Yuwana and Seborg (1982).
The damping ratio ξ is estimated and deduced for a system with τ z = 0 in Söderström and Stoica (1989) Ch. 3 Eq. (3.3g), e.g., where the relative overshoot v is estimated as A proof of the expression of the damping ratio in Eq.
(9) is given in App. F. The steady-state gain of the closed loop system is defined as By putting the derivative of Eq. (5) w.r.t. time ( dy dt = 0) equal to zero, thereafter solving w.r.t. t = t p , we have that For the 1st under-/overshoot we have that Using Eqs. (12) and (14) we obtain Solving Eq. (15) for τ gives Evaluating Eq. (5) at t = t p gives then by using Kr = y ∞ , and solving Eq.
where we have defined that Where we in Eqs. (17) and (18) have used that sin(ωt p + φ) = 1 − ξ 2 , as shown in App. F. Last, we solve Eq. (6) w.r.t. τ z and get that The algorithm is presented and summed up in Alg. 2.1.

Algorithm 2.1 (The JSDR algorithm)
Step 1. From the step response time series data Y and R in Eq. (4) and illustrated in Figs. 2 and 3, calculate the six parameters (t p , y p , r, t u , y u , y ∞ ). An algorithm for step 1 is presented in App. A.
Step 3. Back calculate to the proposed open loop model structure in Eq. (1).
These adjustments propose the improved algorithm. The JSDR algorithm m-file implementation can be found in App. B.
In Jahanshahi and Skogestad (2015) Eq. (A.17) will not work for a case of inverse response, i.e. τ z < 0. This is solved by using Eq. (20) instead, however the drawback is that we need to know a priori if τ z < 0.
Interestingly, an error was found in Eq. (A-5) in Yuwana and Seborg (1982), which estimates the time instant, t p , of the first peak. This error has propagated into wrong estimation of φ in Eq. (A.12), in Jahanshahi and Skogestad (2015). Notice here that Eq. (13) for φ should be used instead of the wrong Eq. (A.12) in Jahanshahi and Skogestad (2015). Note also that this error has no direct impact on the Yuwana-Seborg algorithm, but it impacts the estimate of τ z in Eq. (20).

Proposed alternative algorithm (DR)
An alternative closed loop time transient response to Eq. (5) can delicately be expressed as where ω is the angular frequency, as defined in Eq. (7), and Now, by simply solving Eq. (21) w.r.t. c at t = t p , we obtain that Notice that the above parameter c is related to the phase angle φ (Jahanshahi and Skogestad (2015)) as and as in Eqs. (8) and (13). By solving Eq. (22) w.r.t. τ z , we have that Eqs. (25) or (26) is more preferable than Eq. (20), partly because Eq. (20) is observed to produce complex solutions, hence, τ z = real(τ z ) is added in the JSDR algorithm. DR will not have the problem of having to know if τ z < 0, in prior, which gives DR an edge over JSDR.
The DR algorithm m-file implementation can be found in App. C. Note that the difference from JSDR in Alg. 2.1 is the estimating of τ z .
Notice that an alternative formula for the time instant of the first extremum Eq. (12) is found by solving the derivative of Eq. (21) w.r.t. time. This gives The time instant of the first under-/overshoot is found to be This gives the same as Eq. (15),which is obtained by putting Eq. (27) into Eq. (28) and solving for t u . Solving Eq. (27) w.r.t. c we obtain We also propose an alternative DR2 algorithm, using c from Eq.(29) instead of Eq. (23). The DR2 algorithm m-file implementation can be found in App. D. Notice however, that this alg. DR2 is more sensitive to noise than both JSDR and DR.

The propagation of errors based in the Yuwana-Seborg article
Notice that in Yuwana and Seborg (1982) a first order plus time delay model, i.e.
as the assumed open loop model is considered. In Yuwana and Seborg (1982), Eqs. (9) and (10) are both wrong. These equations can also be found in Hapoglu et al. (1998).
However, the correct equations are and where σ 1 = ξ √ K + 1 + ξ 2 (K + 1) + K − 1 and It might be on purpose that "-1" was removed in Eq. (33) for avoiding complex solutions. This error/approximation was also briefly noted in Jutan and Rodriguez (1984) and Taiwo (1993). • JSDR and DR are equipped with the following function, y = f iltf ilt(b f , a f , x), with motivations for the incoming example. Here, 'filtfilt' is a function for zero-phase digital filtering, and b f and a f are from the function 'butter'. Both these functions are contained in the Signal Processing Toolbox in MATLAB. This filter is however manually tuned for the incoming example.

Further details and discussions
• An additional subspace algorithm, denoted DSR ry, is presented under App. E. This algorithm is using the DSR algorithm on the closed loop data R and Y, then back calculating to the assumed open loop model.
• Note that the presented algorithms may automatize various existing tuning methods, e.g. the 'Good Gain method' presented in Haugen (2012). Eqs. (1) and (35) may be formulated as a discretetime closed loop system,

Numerical examples
where the process noise v k and measurement noise, w k are chosen as white noise with covariance matrices E(v k v T k ) = 0.5 2 I 2 and E(w k w T k ) = 0.05 2 I 1 , respectively.
Two different reference signals, r k , are used, i.e. PRBS for the subspace algorithms/PEM and a simple step for the transient algorithms.
The For comparing the results from the Monte Carlo simulations, we consider the following performance criterion, which is given by Eq. (B.67b) in Söderström and Stoica (1989): is the true parameters, M simulations and N samples.
Results from the Monte Carlo simulations are shown in Figs. 5, 6 and Tab. 1. Our proposed algorithm, DR, is far more robust than the JSDR algorithm. Furthermore it shows that the subspace algorithms and PEM render themselves superior over the transient algorithms. PEM has the best performance while DSR ry is the runner-up candidate just before DSR e. Note that DR2 alg. is dropped in this example, reason being rather unstable results.

Example 3.2 (van de Vusse chemical reactor)
A chemical isothermal reactor (van de Vusse (1962)) is studied in this example. Consider the following closed where the reaction rate coefficients are k 1 = 50, k 2 = 100 and k 3 = 10. u is the feed flow rate, y is the concentration of the product out the outlet. The concentration of the byproduct into the reactor, v, is treated as an unknown constant or slowly varying disturbance with nominal. value v s = 10.
A sufficient large proportional gain is chosen, K p = 140, i.e. large enough for stabilizing the closed loop system and giving a relativity large overshoot preferable for the transient algorithms.
Choosing a reference signal r = 1.02 gives the steady state x s = [2.1046, 0.8847] T . The closed loop system is then simulated on time 0 < t < 0.20 (Fig. 7) with a sampling time equal 0.0001, given the initial state x(t = 0) = x s and a reference step for the transient algs. and a PRBS for the subspace algs.
However, it should be noted that the generated data is detrended, i.e. Y := Y − Jy s , R := R − Jr and U := U − Ju s , where y s = x s 2 , J is a vector of ones of correct length and N is the length of e.g. Y .
In Tab. 2, the DR2 algorithm estimates a model which is closest to the Jacobian estimate, based on y s , u s , x s .
We compare the identified models, generated from the following algs., JSDR, DR, DR2, DSR e and DSR ry, over M = 100 different validation sets, i.e. M different PRBS reference realizations.
The performance is measured by the criterion, whereŷ i is the true output and y i is the simulated output from the identified models. Results from the simulations are shown in Tab. 3. In this particular example it is shown that the JSDR alg. is the best performing of the transient algs., however the subspace algorithms gives the best models.

Discussion
• Notice that by defining v as in Eq. (10) and evaluating the time response, Eq. (5), it can be shown that the expression for ξ becomes as in Eq. (9). This was however also pointed out in the appendix in Yuwana and Seborg (1982).
• Notice that the back calculation step to developing the open loop model is not unique, as also described in Secs. 2.2 and 2.3.
• Notice, relative small sampling time dt have to be used for obtaining accurate estimates.
• As demonstrated in Ex. 3.1 the transient identification methods need a large number of observations compared to the classical subspace and prediction error methods. Even of this, the quality of the estimation results is poor compared to PEM/DSR as illustrated in Tab. 1.

Concluding Remarks
A couple of errors in Yuwana and Seborg (1982) have been documented, one of which caused implications for the algorithm presented in Jahanshahi and Skogestad (2015).
A new alternative formulation for the algorithm presented in Jahanshahi and Skogestad (2015) is proposed and compared with Monte Carlo simulations. The new algorithm DR is found to outperform the JSDR algorithm, however these transient response algorithms are rather inferior in comparison to the highly respected subspace algorithms, i.e. DSR/DSR ry (Di Ruscio (1996)), DSR e (Di Ruscio (2009)) and the classical prediction error method (PEM), e.g. Ljung (1998).
A. Function for estimating the 6 variables in Figure 2. The peak at the second extremum t = t p +t u is given by where we have used that sin(φ + π) = − sin(φ) for an angle φ. From this we evaluate the numerator and denominator in Eq. (10) as Substituting Eqs (49) and (50) Hence, the relative overshoot v is independent of the term D (Eq. (6)) and hence independent of the system zero t z . Using Eq. (15) for t u we obtain From Eq. (52) we have the equation Solving Eq. (53) for ξ 2 we obtain ξ 2 = ln 2 (v) π 2 + ln 2 (v) .
Assuming 0 ≤ v ≤ 1 we have to use the positive root and hence where 0 ≤ v ≤ 1 is the relative damping.