Next Article in Journal
Curvature Invariants of Statistical Submanifolds in Kenmotsu Statistical Manifolds of Constant ϕ-Sectional Curvature
Next Article in Special Issue
From Identity to Uniqueness: The Emergence of Increasingly Higher Levels of Hierarchy in the Process of the Matter Evolution
Previous Article in Journal
Category Structure and Categorical Perception Jointly Explained by Similarity-Based Information Theory
Previous Article in Special Issue
Relating Vertex and Global Graph Entropy in Randomly Generated Graphs
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Free Final Time Input Design Problem for Robust Entropy-Like System Parameter Estimation

Faculty of Computer Science, Bialystok University of Technology, Wiejska 45A, 15-351 Bialystok, Poland
Entropy 2018, 20(7), 528; https://doi.org/10.3390/e20070528
Submission received: 3 June 2018 / Revised: 6 July 2018 / Accepted: 12 July 2018 / Published: 14 July 2018
(This article belongs to the Special Issue Entropy: From Physics to Information Sciences and Geometry)

Abstract

:
In this paper, a novel method is proposed to design a free final time input signal, which is then used in the robust system identification process. The solution of the constrained optimal input design problem is based on the minimization of an extra state variable representing the free final time scaling factor, formulated in the Bolza functional form, subject to the D-efficiency constraint as well as the input energy constraint. The objective function used for the model of the system identification provides robustness regarding the outlying data and was constructed using the so-called Entropy-like estimator. The perturbation time interval has a significant impact on the cost of the real-life system identification experiment. The contribution of this work is to examine the economic aspects between the imposed constraints on the input signal design, and the experiment duration while undertaking an identification experiment in the real operating conditions. The methodology is applicable to the general class of systems and was supported by numerical examples. Illustrative examples of the Least Squares, and the Entropy-Like estimators for the system parameter data validation where measurements include additive white noise are compared using ellipsoidal confidence regions.

1. Introduction

System identification is typically carried out by perturbing processes or plants under operation and use experimental data to construct the model of the dynamic system. The main objective is to find a set of dynamic models that describe important properties of the true system [1]. The fundamental task in system identification is to excite the system of interest using an informative input and build the model of the system with maximum pertinence [2,3]. The problem of the optimal input signal design is typically solved by minimizing an a priori selected norm of the Fisher information matrix with respect to an appropriate experimental setup [4]. The prediction error methods (PEM) are a wide collection of the parameter estimation methods that minimize a weighted norm of the prediction error [2]. The identification experiment can be executed in both closed and open loop conditions and could be utilized for arbitrary model parameterizations. Improper experiment conditions can cause performance degradation of the control loop. It has been reported that about 80% of the designed control loops do not guarantee the acceptable performance assessment [5].
That is why some authors introduced the idea of the performance degradation minimization instead of the variance minimization of the estimated parameters. The robust control identification considers the uncertainty of the estimated model on the designed closed-loop system performance [6]. System identification for robust control allows comparing the performance of the unknown real system loop with a controller tuned using an identified plant model. The least-costly identification experiment for control, where the main goal is to design an experiment that ensures a small enough uncertainty region but still provides an acceptable performance of control, was proposed in [7,8]. It was found that, during the advanced control loop designs, model construction absorbs approximately 75% of the costs [9].
The plant-friendly input design is classified as the application-oriented methodology. The aim of such an identification experiment is to find a trade-off between the minimal disruption to the normal operation of the system, and the most precise identification experiment [10,11]. There have been some reports that plant friendliness constraints often disturb a precise model parameters estimation while a set of harmonically related sinusoids with high peak-to-peak values can destroy an identified model [3,12]. For this reason, safer excitation signals which ensure more precise model identification are presented in [13,14]. The issue of optimal input design in the economic framework, where the cost based on plant friendly constraint is minimized, was developed in [15]. One of the current trends in accordance with application-oriented input design is the use of the model predictive control (MPC) technique [16,17]. This formulation is based on the input design procedure to obtain an acceptable control performance that still provides revealing data for system identification [18]. The main idea is to choose the spectral density of an input signal that ensures that estimated parameters are acceptable while the experimental cost is minimized. For an overview of application-oriented input design in application to system identification, see the survey [19]. In the exact model identification for control purposes, the optimal experiment design methods should be considered. For this reason, Entropy-based optimal experiment design method for discrimination between competing models was presented in [20]. The approach of model discrimination was based on the expected Shannon entropy reduction of the Bayesian model weights uncertainty.
Considering the automatic control aims and objectives for model identification, the outlying data have a critical impact on model parameters to be estimated. The most utilized prediction error estimator methods are Least Squares (LS) and Weighted Least Squares (WLS) where the sum of the squared residuals is minimized [21]. Different robust formulations [22] where the penalty function is minimized with respect to the overall distribution of residuals are Least Median of Squares (LMS), Least Trimmed Squares (LTS), and Reweighted Least Squares (RLS). The Maximum Likelihood (ML), Minimum Entropy (ME), and Generalized Maximum Entropy (GME) approaches for robust parameter identification, which guarantee robustness subject to regression models, were presented in [23]. The novel prediction error parameter estimation method named Least Entropy-Like (LEL) estimator was developed in [24]. This algorithm is based on correctly established penalty function and was built according to the Gibbs entropy definition. In a previous paper [25], the spectrum approximation issue based on the idea of optimal prediction was presented. The THREE-like approach minimizes a divergence factor subject to a given spectral density. The contribution of this work was to define a new divergence family which compares two spectral densities with respect to an optimal prediction task. In this approach, the output covariance of a bank of filters was utilized to obtain information on the input spectrum power. The interpretation of the dual problem of the THREE-like approach is shown as a new parametric spectral approximation problem wherein the optimal result of the spectral density is closest to the correlogram. It has been shown that two particular THREE-like solutions are equivalent to the prediction error identification method [26]. In the above paper, an important connection between time and spectral domain entropy rates was presented in details.
After a brief review of current literature, there is still an interest in parameter identification from real-life data [27]. The methodology presented in this paper can be successfully used in continuous and discrete-time state-space model identification including linear and nonlinear models with higher orders. In a previous paper [28], a first-order linear time-invariant system case study of plant friendly input signal design with respect to the D-efficiency was presented. The economical and plant friendly input design for linear discrete-time system identification, where the goal was to minimize the cost incurred during the experiment with respect to plant friendly constraints, can be found in [15]. Numerical results for second-order linear torsional spring system identification were introduced in [29]. The optimal input design problem for parameter estimation in a nonlinear water tank system through a model output sensitivity minimization subject to its parameters was solved in [30]. However, all the existing works in the literature consider the optimization methods of designing the excitation signals with fixed final time conditions only. In contrast, the present article discusses the formulation and the solution scheme of a free final time optimal input design problem. For this purpose, the state-space equation is augmented by an extra state variable representing free terminal time scaling factor (i.e., according to [31]). The contribution of this study is to use the free terminal time inputs to examine the economic aspects between the imposed constraints on the input signal shape, and the parameter estimates while undertaking a robust identification experiment utilizing an implementation of the LEL algorithm. The constraints imposed on input signal solution should allow one to attain a slight information loss from the plant whose operating point is disturbed in the safest way.
This paper is structured as follows. In Section 2, the problem statement of the constraint optimal input design for parameter estimation of a system is described. The Least Entropy-Like estimator for the robust system identification is presented in Section 3. In Section 4, D-optimal input signal design and the transcription for free final time problem is derived. The problem reformulation for free final time constraint input design is shown in Section 5. The results of simulation experiments for a linear time-invariant model case study are presented in Section 6. Finally, concluding comments are made in Section 7.

2. Problem Statement

Consider below system described as:
y ( t ) = f ( u , t , θ , v ) ,  
where y(t), u(t), v(t) and θ are, respectively: output, input, noise, and parameter of the system. System identification is the process of building an accurate dynamic mathematical model of a system from experimental data and a priori plant knowledge. The precision of the model parameter estimates depends on the selection of an optimal perturbation signal. In the design of optimal input signals for a plant model parameters calculation, an appropriate scalar norm of the Fisher information matrix (FIM) must be chosen as the objective criterion. The FIM is defined as follows:
M = E [ θ ln p ( y | θ ) ( θ ln p ( y | θ ) ) T ] .  
The D-optimality norm for which determinant of the Fisher information matrix detM is maximized or det(M−1) is minimized is often used. Chosen measures of optimal design performance could be found in [32]:
  • A-optimality (tr(M−1)) minimizes the total variance of the parameters estimates.
  • E-optimality (λmax(M−1)) minimizes the variance of the maximum eigenvalue of M−1.
  • D-optimality minimizes the generalized variance of the parameters and minimizes the volume of the ellipsoidal confidence region of parameter estimates with respect to the input.
To obtain optimal input signal, an unbiased estimator of θ should be considered. In this case, the covariance of the parameter estimates is given by the Cramer–Rao inequality, viz., the inverse of the FIM. Then, the covariance of the estimate θ ^ can be defined as:
cov ( θ ^ ) = E [ ( θ ^ θ ) ( θ ^ θ ) T ] M 1  
The single weighted cost function method for time domain input signal design was presented in a previous paper by the author [28]. While the choice of the experiment norm is essential, the inputs designed based on some performance criteria may not be suitable for the plant excitation [2]. The input signal used for a plant model excitation should, at the same time, fulfill two requirements: the acceptable precision of the model parameter estimates and the system should be perturbed in the safest way. These conditions can be met using an approach, which is based on the notion of the D-efficiency [28]. The D-efficiency is often expressed as a percentage scale and may be considered as a measure of the sub-optimality of any synthesized input signal:
E D ( e ) = 100 % { det ( M ( e ) ) det ( M ( e * ) ) } 1 k ,  
where k denotes the number of parameters to be identified, and e* indicates the D-optimal design. Regarding the analysis presented in t [28], we impose the inequality scalar constraint on the D-efficiency formulation. In general, the optimal input design task is formulated through maximization of the FIM determinant and should take into account the equality and inequality constraints imposed on the conditions. The experiment performed in such a way should yield safe behavior of the perturbation signal. Some of these constraints could be defined as:
  • D-efficiency inequality constraint can be expressed according to (4).
  • Amplitude constraints on inputs, outputs or state variables in a form:
    u min u ( t ) u max .  
  • Energy constraints imposed on inputs:
    0 T u T ( t ) u ( t ) d t E .  
The goal of this study was to design the constraint input signal with free final time conditions, which is then used to the robust identification experiment. In this case, the performance index is formulated through the minimization of the free terminal time scaling factor, subject to the D-efficiency constraints as well as input energy constraints. In that way (i.e., by taking into account such a set of constraints to optimal input design), we can obtain the suboptimal input signal, which is safer for system identification purposes.

3. System Identification Method

In system identification for control purposes, outlying data have a critical impact on plant model development including precise parameter estimation. In this study, the least squares and the novel Entropy-Like estimators were compared for plant model parameters identification purposes. For these reasons, the basic concept of the prediction error method for system identification is presented in this section. The algorithm is based on the specific objective function minimization with respect to prediction error residuals. The concept of the cost function formulation was motivated by the Gibbs entropy definition. A more detailed overview of the Entropy-Like method (LEL) was presented in [24].
To solve the optimization task presented above, the following system was considered:
y i = f ( x 1 , x 2 , , x l , θ r ) + ε i , i = 1 , 2 , , l ,  
where θr is a parameter vector to be estimated, yi is the output samples sequence, and εi is the Gaussian white noise with restricted variance. The model prediction error estimators are obtained utilizing the regression residuals:
r i = y i y ^ i ,  
where y ^ i are the approximated outputs, i.e., y ^ i = f ( x 1 , x 2 , , x l , θ ^ ) . The Least Squares method (LS) for regression analysis is standard and very popular approach. The performance index minimizes the sum of squared residuals of the fit (i.e., prediction error estimator) is:
θ ^ L S = Arg min θ i = 1 N r i 2 = Arg min θ ( r T r ) ,  
where r is the residual vector of the fit r = f ( r 1 , r 2 , , r N ) . A different algorithm, which is robust for model parameter identification, is based on penalty function minimization.
The Least Entropy-Like estimator (LEL) was developed utilizing the concept of Gibbs entropy. The main objective of this method is to analyze the global dispersion measure of the residuals fit. The prediction error estimator built according to Equation (8) is as follows:
S = j = 1 N r j 2 ,  
the comparative squared residuals can be presented as:
if     S 0 ​​​    then     s i = r i 2 j = 1 N r j 2     where     s i [ 0 , 1 ] , i = 1 N s i = 1 ,  
according to the reasoning presented in [24], the cost function Φ based on normalized entropy was chosen as:
Φ = { 0                                             for     S = 0 1 log N i = 1 N s i log s i     for     S 0 .  
Then, the Least Entropy-Like (LEL) estimator defined as a measure of dispersion of the relative squared rest values is given by:
θ ^ L E L = Arg min θ Φ ,  
Entropy-like formulation (12) affects the values of the unknown parameters θ through the predictive error residuals. The LEL estimator (13) is robust with respect to outliers because the objective function is minimized subject to relative squared errors variability. It could be noticed that the cost function Φ is nonlinear and may not provide the unique minimum subject to the unknown parameters θ. Considering the basic algorithm properties, one should initially examine the Least Squares quality fit. As presented in [24], the LEL estimator can be numerically computed from an initial condition value close to the real parameter value. In the experimental part of this paper, the LS and the novel LEL estimators would be compared for plant model parameters identification purposes.

4. Optimal Input Design Problem

In the paper, the problem of synthesizing the constraint optimal input with free terminal conditions for system identification is considered. The general idea is to define a nominal period [0, Tf] and to replace free final time with fixed final time problem utilizing a scaling factor as an augmented state variable, which scales the duration of the time interval. We solve this problem using the transcription of the below optimal control formulation into a similar optimal control task represented in the Lagrange form with the set of constraints. To verify the suitability of this technique to the model parameter identification, a first-order time-invariant inertial system is considered:
x ˙ ( t ) = a x ( t ) + b u ( t ) ,    x ( 0 ) = x 0 , y ( t ) = x ( t ) + v ( t ) ,  
where x(t) is a state variable, u(t) is a control function, y(t) is a measurement, a and b are constant model parameters and v(t) is a zero mean Gaussian white noise process as follows:
E [ v ( t ) ] = 0 , E [ v ( t ) v T ( τ ) ] = R δ ( t τ ) = σ n 2 δ ( t τ ) .  
The fundamental principle of system parameter estimation is to maximize the sensitivity of the state variable to the unidentified parameters [1]. The motivation for such an approach is the Cramer–Rao definition, which gives a lower bound for the variance of an unbiased parameter to be estimated. Applying the above definition to input design purposes, we calculate the parameter estimate which has a tendency of getting lower for optimal input:
cov ( [ a , b ] ) M 1 .  
The FIM for the inertial state-space model (14) can be formulated as follows:
M ( T ) = 0 T X θ T R 1 X θ d t = 1 σ n 2 0 T [ x a x b ] [ x a x b ] d t = 1 σ n 2 0 T [ x a 2 x a x b x a x b x b 2 ] d t ,  
where xa =x/a, xb =x/b, and R is 2 × 2 matrix given by:
R 1 = [ R 1 0 0 R 1 ] = 1 σ n 2 [ 1 0 0 1 ] .  
Substituting Equation (17) into Equation (16) yields:
cov ( [ a , b ] ) σ n 2 M ( T ) .  
In the considerations which follow, it was assumed that σn = 1 to obtain an optimal input signal for system parameters identification where measurements do not include additive white noise. To maximize the FIM determinant, let us define the augmented state vector given by [1]:
x ˙ a = x + a x a ,     x a ( 0 ) = 0 ,  
x ˙ b = a x b + u ,     x b ( 0 ) = 0 .
The Fisher information matrix to the inertial model (14) has the following form:
M ( t ) = [ m 11 ( t ) m 12 ( t ) m 21 ( t ) m 22 ( t ) ] .  
To design an optimal input signal with free terminal time conditions, the optimal control problem solver RIOTS_95 was utilized [33]. The Matlab toolbox Riots allows solving a very large class of finite-time optimal control problems that includes: trajectory and end-point constraints, variable initial conditions, free final time tasks and problems with cost functions endpoint. The objective function to be minimized can be formulated in Bolza form as:
J = g ( x ( t 0 ) , x ( t f ) ) + t 0 t f l ( x , u , t ) d t ,  
in respect to the plant dynamics, and with the initial condition:
x ˙ ( t ) = h ( x , u , t ) , ​​   x ( t 0 ) = x t 0 , t [ t 0 , t f ] ,  
subject to the constraints:
u ( t ) u min ( t ) , u max ( t ) , t [ t 0 , t f ] ,  
x ( t 0 ) x min ( t 0 ) , x max ( t 0 ) ,  
l t c ς ( t , x ( t ) , u ( t ) ) 0 , ς q tc , t [ t 0 , t f ] ,  
g e i c ς ( x ( t 0 ) , x ( t f ) ) 0 , ς q eic ,  
g e e c ς ( x ( t 0 ) , x ( t f ) ) = 0 , ς q eec .  
where x is the state-space vector, t ∈ [t0, tf] denotes time duration, q = {1, …, q} and l, g, and h are a priori linear or nonlinear functions. The functions g(·,·) and l(·,·,·) with indexes tc, eec, and eic are trajectory constraint, endpoint equality constraint and endpoint inequality constraint, respectively.
The free terminal time problem can be included in the form of an optimal control problem by augmenting the state equations by further state variables, i.e., one extra state for each independent problem. According to the reasoning presented in Section 2 of the user manual [33], the general idea is to define a nominal time interval [0, Tf] and to replace the free final time problem with the fixed final time case utilizing the free final time scaling factor as an augmented state variable, which scales the duration of the time interval. For this reason, the scale factor and the scaled time are expressed by extra states which enables minimization over initial value of the further states to fit the scaling.
Assuming that the state space differential equation is described by:
x ˙ = h ( t , x , u ) ,  
the cost function has the following form:
J = g ( x ( T ) ) + 0 T l ( t , x , u ) d t ,  
where x(t) is the state space vector, u(t) is the state input vector, and T denotes the free terminal time of the experiment.
Including two extra state variables xn+1 and xn+2 to (30), the free terminal time problem can be modified into the similar fixed final time optimal control problem with an augmented state vector:
( x ˙ ( t ) x ˙ n + 1 x ˙ n + 2 ) = ( x n + 2 h ( x n + 1 , x , u ) x n + 2 0 ) ,  
where t ∈ [0, Tf] and the objective function can be written as:
J = g ( x ( x n + 2 T f ) ) + 0 T f l ( x n + 1 , x , u ) d t ,  
where t ∈ [0, Tf], xn+2 is the duration scale coefficient to be minimized, xn+1 = txn+2 denotes free termination time, and Tf is the fixed termination time chosen arbitrarily. When considering the autonomous dynamic systems, the extra state variable xn+1 is not obligatory. Therefore, the autonomous free final time problem can be solved by augmenting state equations considering only one state variable representing the free final time scaling factor.

5. Problem Reformulation

As described in [29] with reference to the single weighted cost function method, the derivations for inertial system dynamics are given below. The augmented state Equations (20) and (21), and the FIM matrix elements containing an extra state variable representing free final time scaling factor ζ are as follows:
x 1 = x ; x ˙ 1 = x 7 ( a x 1 + b u ) ; x 1 ( 0 ) = x 10 ; x 2 = x a ; x ˙ 2 = x 7 ( x 1 + a x 2 ) ; x 2 ( 0 ) = 0 ; x 3 = x b ; x ˙ 3 = x 7 ( a x 3 + u ) ; x 3 ( 0 ) = 0 ; x 4 = m 11 ; x ˙ 4 = x 7 x 2 2 ; x 4 ( 0 ) = 0 ; x 5 = m 12 = m 21 ; x ˙ 5 = x 7 x 2 x 3 ; x 5 ( 0 ) = 0 ; x 6 = m 22 ; x ˙ 6 = x 7 x 3 2 ; x 6 ( 0 ) = 0 ; x 7 = ζ ; x ˙ 7 = 0 ; x 7 ( 0 ) = ζ 70 .
The proposed method can be applied to a general class of systems given by the state-space representation of any order. Finally, the optimal control problem based on Bolza canonical formulation, which minimizes the objective function (23) with respect to D-efficiency equality constraint (4) as well as input (25) and trajectory (27) constraints, is:
J ( u ) = ζ + q 0 T u ( t ) T u ( t ) d t ,  
with respect to:
1 u ( t ) 1 ,    t [ 0 , T ] , 0.1 ζ 10 , x 1 ( t ) 1 ,     t [ 0 , T ] , ( x 4 ( T ) x 6 ( T ) + x 5 2 ( T ) ) = D ,  
where D is the desired D-efficiency constant, q is an input energy factor and ζ70 is an initial state condition to be optimized. It should be noted that an additional constraint was imposed on the state variable x1(t) to enable unexpected changes of the control signal u(t), which is restricted to the interval [−1, +1]. Using the proposed methodology defined by Equations (34)–(36), the optimal solution for free final time is tf = .

6. Numerical Results

To solve the issue presented above, the RIOTS_95 toolbox [33] for solving optimal control problems can be adapted. This software is implemented in Matlab as a separate module and has tools for solving constrained optimal control problems with fixed or free final time conditions.
Constrained optimal inputs for the first-order, linear time-invariant (LTI) model of the system were then computed for the assumed initial values of parameters: a = −1, b = 1, and nominal time duration t = [0, 10] s, using sequential quadratic programming (SQP) algorithm. The initial state conditions of the inertial model were selected to be x1(0) = 5, x7(0) = 1, and the initial value of the input signal was set as u(0) = 1. The free final time scaling factor ζ which scales the duration of the time is optimized from the interval 0.1 ≤ ζ ≤ 10, so the time duration could be varied from 1 to 100 s. The numerical results were computed using the fixed step-size fourth-order Runge–Kutta integration method with grid period of 0.2 s. The expression for the cost function, given by Equation (35), can be presented as:
J ( u ) = J 1 + q J 2 ,  
where J1 denotes the free final time scaling factor ζ, and J2 is the integral of the squared input signal.

6.1. Free-Final Time Constraint Input Design

The D-optimal input signal received when there was no constraint on the input energy value (i.e., the coefficient q ≈ 0 in the equality (35) is displayed in Figure 1a. It corresponds to the D-optimal experiment e, where the desired value of the FIM determinant is obtained (according to (4)) in such a manner that Deff = 90% of its optimal value. The suboptimal input signals obtained for different desired values of the input energy factor q and D-efficiency constant value Deff = 90% are displayed in Figure 1c,d.
The D-optimal perturbation signal received when there was no constraint on the input energy component (i.e., for J1 = 0.88, qJ2 = 1.0 × 10−4 and tf = 8.79 [s]) is shown in Figure 1a. The input energy factor was increased (Figure 1b) to obtain the suboptimal input signal, which corresponds to performance index components values: J1 = 0.92, qJ2 = 5.00 and tf = 9.23 [s]. For comparison, Figure 1c shows the suboptimal input signal, which correlates with the objective function values: J1 = 0.97, qJ2 = 8.70 at the final time level of tf = 9.67 [s]. Figure 1d contains the graphical display of the suboptimal excitation received for the cost function integrants values: J1 = 1.01, qJ2 = 33.44, where time duration was tf = 10.11 [s].
As we can see, when the desired value of input energy factor increases, the shape of the optimal excitation considerably changes. While for the optimal experiment (in the sense of Equation (4)) there are the rapid changes of the input, the control signals obtained for Deff < 100% are safer for system identification purposes until the FIM determinant is not dominated by the input energy component of the minimized performance index. The comparison of the performance index components obtained for increasing values of the input energy factor and for decreasing values of D-efficiency from the interval [100%, 80%] of its maximum value are presented in Table 1, Table 2 and Table 3.
As could be noticed based on the presented method (see Table 1, Table 2 and Table 3), when the desired value of the input energy factor increases, the D-optimal signal duration also increases. When the required value of the D-efficiency from the interval [100%, 80%] decreases, considered signals duration also decreases. As we can see (Table 1), the optimal input signal duration for inertial system identification (i.e., J1 = 1, q = 1 ×10−6 and Deff = 100%) is equal to 10 s.

6.2. LS and LEL Estimators for LTI Model Identification

The D-optimal input signals u(tf), computed as solutions of the free final time optimization problems (35) and (36), were then utilized as excitations in the plant model parameter estimation procedure. The physical system (14), used in system identification procedure, can be described by the following single input–single output state space model:
x ˙ 1 ( t f ) = a x 1 ( t f ) + b u ( t f ) + v ( t f ) ,    x 1 ( 0 ) = x 10 , y ( t f ) = x 1 ( t f ) ,  
The scheme shown in Figure 2 presents the process of the plant model parameter estimation: we excite the system input using u(tf), and we collect data on its output y(tf).
A disturbance signal with different variance from the interval 0.0 ≤ σ2 ≤ 0.7 is added to the control input to the system. The model of the plant (38) depends on a vector of unknown parameters θ = [a, b]T and the aim of such an experiment is to estimate unknown model parameters values which should be the most similar to the true values of the plant parameters. The difference between the output of the plant y(tf) and the output of the model ym(tf) was minimized. The initial state condition of the inertial model was selected from the interval −5 ≤ x1(0) ≤ 5 and the experiment duration depends on chosen D-optimal signal according to Table 2. Numerical results were obtained utilizing the Nelder–Mead simplex method.
The distribution of the model parameters a and b obtained as the results of the optimization tasks using conventional LS and robust LEL estimators with different control inputs are shown in Figure 3.
Eighty-eight computations were done when the plant model starts from various initial state condition, and the additional noise disturbing the system input has a different variance. Figure 3a contains the ellipsoidal confidence region (black dotted line) with the input signal obtained for the minimal value of input energy factor (i.e., q ≈ 0 and Deff = 90%), while Figure 3b shows the results (for the same values of initial states and noise variance) with excitation signals that were computed when an input energy coefficient increases its value and was selected as q = 0.10 (green dashed line) and 0.40 (red dash-dot line). To compare results, Figure 3 contains the graphical display of the ellipsoidal confidence region of parameter estimates, where the system (Figure 2) was perturbed utilizing a step input signal (blue solid line). The comparison of the ellipsoidal confidence regions of the plant model parameter estimates indicates some similarities. The D-optimal input signal, calculated for q ≈ 0, causes the minimal time duration of the parameter identification experiment and a minimal volume of the ellipsoidal confidence region of parameter estimates. When the value of the input energy factor increases, the area under the curve increases its size for the same initial conditions and noise variance values. Increasing the desired ratio of input energy constraint yields the increase in the input signal duration, but the excitation is safer for the plant. In such a manner, we avoid rapid switching of the excitation signal in the real identification experiments. The prediction error LS, and the relative squared error LEL estimators were compared based on maximum, average, and minimum residuals of the parameter values (Table 4 and Table 5).
As could be noticed based on Table 4 and Table 5, a properly implemented penalty function of the LEL algorithm yields, in most cases, more precise estimates for different values of the energy constraint, D-efficiency constraint, and the experiment duration. The comparison of the estimates for parameters a and b is shown in Figure 4.
The function Φ defined by Equation (12) depends on parameters a and b by means of relative squared residuals. The main idea of the LEL estimator (13) is to make most of the residuals striving for zero value or to make the relative squared residuals equally distributed according to minimization criteria. Data points related to large residuals are called the outlying data points. The proposed method is robust with respect to outliers because the cost function to be minimized is a measure of the relative squared errors variability. In the above simulation experiment, the outlier data occurred for inertial system initial condition x1(0) = 1. The reason for this obstruction is the excitation signal (i.e., unit step-like input) which is not able to unbalance the inertial system with respect to the initial condition equal to 1. Therefore, data outlier points were removed from the set of the parameter estimates. Finally, it should be noticed that there is no guarantee for relative squared residuals formulation to have a unique solution with respect to the parameters. Thus, the minimization should be performed very carefully with special attention given to the initialization of the parameters.

7. Discussion

The novelty of this work was to design the free final time input signals, with constraints on input move size and D-efficiency value which are then used in the system identification experiments. The objective of such an experiment was to minimize the extra state variable representing the free final time scaling factor, subject to a carefully selected set of constraints. Some of these signals obtained as the solution of the free final time optimization task were then used as inputs in the plant model parameter estimation procedure. The parameter identification experiments were performed in the presence of a white noise affecting the system input. The imposed constraints, which quantify the plant trajectories, allow one to obtain safer excitation signals while still providing acceptable confidence regions of the plant model parameter estimates.
Another objective of this work was to find the relationships between the imposed constraints on the input signal design and the experiment duration in the real-time parameter identification experiments. The results of the above experiments confirm our assumption that the input signal obtained for input energy factor value q ≈ 0, yields the minimal duration of the identification experiment and a minimal volume of the ellipsoidal confidence region of parameter estimates. When the D-efficiency constraint percentage ratio decreases, then the experiment duration is reduced as well as the economic cost of identification under normal operation. However, safer input signals obtained for different required values of the input energy factor yield the relative accuracy degradation of the plant model parameter estimates (observed as an increased volume of the confidence ellipsoid regions).
An identification cost for industrial purposes is difficult to obtain because it must contain actual costs during the regular production process. The signals used to system perturbation or resulting outputs of the identified system can affect other production processes. In the presented method, the cost is quantified in terms of departure from the nominal operational policy. During the closed loop identification experiment, the framework discussed in this paper is utilized to design a reference signal instead of the input signal to the open-loop model. The solution of the identification task appears to be robust with respect to the initial conditions variation and was expressed as a non-convex problem.
It was shown that applying the proposed method to system identification purposes yields an interesting empirical solution. The results obtained in the identification experiments prove that there is a trade-off between the experiment duration (i.e., closely related to input signal duration) and the accuracy of parameter estimates which depends on the friendliness of the input signal to be used. Thus, the cost of the identification experiment should be considered either as the experiment duration or as a measure of the performance deterioration.

Funding

This research was funded by [Politechnika Białostocka] grant number [S/WI/3/18].

Acknowledgments

I am deeply indebted to M. Świercz, Faculty of Electrical Engineering, Bialystok University of Technology, for his guidance and necessary help during the research. The present study was supported by a grant S/WI/3/18 from the Bialystok University of Technology and founded from the resources for research by the Ministry of Science and Higher Education.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Kalaba, R.; Spingarn, K. Control, Identification, and Input Optimization; Plenum Press: New York, NY, USA, 1982; pp. 225–299. ISBN 978-1-4684-7662-0. [Google Scholar]
  2. Ljung, L. System Identification: Theory for the User; Prentice-Hall: Englewood Cliffs, NJ, USA, 1999; pp. 358–406. ISBN 0-13-881640-9. [Google Scholar]
  3. Pintelon, R.; Schoukens, J. System Identification: A Frequency Domain Approach, 2nd ed.; John Wiley & Sons: New York, NY, USA, 2001; pp. 151–224. ISBN 978-0-470-64037-1. [Google Scholar]
  4. Mehra, R. Choice of Input Signals. In Trends and Progress in Systems Identification; Eykhoff, P., Ed.; Pergamon Press: New York, NY, USA, 1981; pp. 305–366. [Google Scholar]
  5. Hugo, A.J. Process Controller Performance Monitoring and Assessment. Control. Arts Inc. 2001. Available online: http://www.controlarts.com/ (accessed on 28 May 2018).
  6. Gevers, M.; Ljung, L. Optimal experiment designs with respect to the intended model application. Automatica 1986, 22, 543–554. [Google Scholar] [CrossRef]
  7. Bombois, X.; Scorletti, G.; Gevers, M.; Van den Hof, P.M.J.; Hildebrand, R. Least costly identification experiment for control. Automatica 2006, 42, 1651–1662. [Google Scholar] [CrossRef]
  8. Bombois, X.; Hjalmarsson, H.; Scorletti, G. Identification for robust deconvolution filtering. Automatica 2010, 46, 577–584. [Google Scholar] [CrossRef]
  9. Hussain, M. Review of the applications of neural networks in chemical process control-simulation and on-line implementation. Artif. Intell. Eng. 1999, 13, 55–68. [Google Scholar] [CrossRef]
  10. Rivera, D.; Lee, H.; Braun, M.; Mittelmann, H. Plant friendly system identification: A challenge for the process industries. In Proceeding of the SYSID 2003, Rotterdam, The Netherlands, 27 August 2003; pp. 917–922. [Google Scholar]
  11. Narasimhan, S.; Rengaswamy, R. Plant friendly input design: Convex relaxation and quality. IEEE Trans. Automat. Control 2011, 56, 1467–1472. [Google Scholar] [CrossRef]
  12. Steenis, R.; Rivera, D. Plant-friendly signal generation for system identification using a modified simultaneous perturbation stochastic approximation (SPSA) methodology. IEEE Trans. Control Syst. Technol. 2011, 19, 1604–1612. [Google Scholar] [CrossRef]
  13. Rafajłowicz, E.; Rafajłowicz, W. A variational approach to optimal input signals for parameter estimation in systems with spatio–temporal dynamics. In Proceedings of the 10th International Workshop in Model-Oriented Design and Analysis, Łagów Lubuski, Poland, 10–14 June 2013; Springer: Heidelberg, Germany, 2013; pp. 219–227. [Google Scholar] [CrossRef]
  14. Rafajłowicz, E.; Rafajłowicz, W. More safe optimal input signals for parameter estimation of linear systems described by ODE. In Proceedings of the 26th Conference on System Modeling and Optimization (CSMO), Klagenfurt, Austria, 9–13 September 2013; Springer: Heidelberg, Germany, 2014; Volume 443, pp. 267–277. [Google Scholar] [CrossRef]
  15. Kumar, A.; Nabil, M.; Narasimhan, S. Economical and plant friendly input design for system identification. In Proceedings of the European Control Conference, Strasbourg, France, 24–27 June 2014; Volume 48, pp. 732–737. [Google Scholar] [CrossRef]
  16. Potters, M.G.; Bombois, X.; Forgione, M.; Modén, P.E.; Lundh, M.; Hjalmarsson, H.; Van den Hof, P.M.J. Optimal experiment design in closed loop with unknown, nonlinear and implicit controllers using stealth identification. In Proceedings of the European Control Conference, Strasbourg, France, 24–27 June 2014; pp. 726–731. [Google Scholar] [CrossRef]
  17. Larsson, C.A.; Annergren, M.; Hjalmarsson, H.; Rojas, C.R.; Bombois, X.; Mesbah, A.; Modén, P.E. Model predictive control with integrated experiment design for output error systems. In Proceedings of the European Control Conference, Zurich, Switzerland, 17–19 July 2013; pp. 3790–3795. [Google Scholar]
  18. Larsson, C.A.; Rojas, C.R.; Bombois, X.; Hjalmarsson, H. Experiment evaluation of model predictive control with excitation (MPC-X) on an industrial depropanizer. J. Process Control 2015, 31, 1–16. [Google Scholar] [CrossRef]
  19. Annergren, M.; Larson, C.A.; Hjalmarsson, H.; Bombois, X.; Wahlberg, B. Application-oriented input design in system identification. Optimal input design for control. IEEE Control Syst. Mag. 2017, 37, 31–56. [Google Scholar] [CrossRef]
  20. Malakara, N.K.; Knuth, K.H. Entropy-Based Search Algorithm for Experimental Design. AIP Conf. Proc. 2011, 1305, 157–164. [Google Scholar]
  21. Rousseeuw, P.J.; Leroy, A.M. Robust Regression and Outlier Detection; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2003; pp. 21–145. ISBN 0271-6356. [Google Scholar]
  22. Huber, P.J. Robust Statistic; John Wiley & Sons, Inc.: New York, NY, USA, 1981; pp. 43–68. ISBN 0-47141805-6. [Google Scholar]
  23. Golan, A.; Judge, G.G.; Miller, D. Maximum Entropy Econometrics: Robust Estimation with Limited Data; John Wiley & Sons Inc.: Chichester, UK, 1996; ISBN 978-0-471-95311-1. [Google Scholar]
  24. Indiveri, G. An entropy-like estimator for robust parameter identification. Entropy 2009, 11, 560–585. [Google Scholar] [CrossRef]
  25. Zorzi, M. Multivariate spectral estimation based on the concept of optimal prediction. IEEE Trans. Autom. Control 2015, 60, 1647–1652. [Google Scholar] [CrossRef]
  26. Zorzi, M. An interpretation of the dual problem of the THREE-like approaches. Automatica 2015, 62, 87–92. [Google Scholar] [CrossRef] [Green Version]
  27. Stojanovic, V.; Filipovic, V. Adaptive input design for identification of output error model with constrained output. Circuits Syst. Signal Process. 2014, 33, 97–113. [Google Scholar] [CrossRef]
  28. Jakowluk, W. Plant friendly input design for parameter estimation in an inertial system with respect to D-efficiency constraints. Entropy 2014, 16, 5822–5837. [Google Scholar] [CrossRef]
  29. Jakowluk, W. Optimal input signal design for a second order dynamic system identification subject to D-efficiency constraints. In Computer Information Systems and Industrial Management; Lecture Notes in Computer Science; Springer: Cham, Switzerland, 2015; pp. 351–362. [Google Scholar]
  30. Chianeh, H.A.; Stigter, J.D.; Keesman, K.J. Optimal input design for parameter estimation in a single and double tank system through direct control of parametric output sensitivities. J. Process Control 2011, 21, 111–118. [Google Scholar] [CrossRef]
  31. Tricaud, C.; Chen, Y. An approximate method for numerically solving fractional order optimal control problems of general form. Comput. Math. Appl. 2010, 59, 1644–1655. [Google Scholar] [CrossRef]
  32. Atkinson, A.; Donev, A.; Tobias, R. Optimum Experimental Design with SAS; Oxford University Press: Oxford, UK, 2007; pp. 135–147. ISBN 9780199296606. [Google Scholar]
  33. Schwartz, A.; Polak, E.; Chen, Y. Riots a Matlab Toolbox for Solving Optimal Control Problems. Version 1.0 for Windows, May 1997. Available online: http://www.schwartz-home.com/RIOTS/ (accessed on 28 May 2018).
Figure 1. Free final time inputs to inertial model perturbation: (a) optimal input signal for q ≈ 0 and Deff = 90%; (b) suboptimal input signal for q = 0.05 and Deff = 90%; (c) suboptimal input signal for q = 0.10 and Deff = 90%; and (d) suboptimal input signal for q = 0.40 and Deff = 90%.
Figure 1. Free final time inputs to inertial model perturbation: (a) optimal input signal for q ≈ 0 and Deff = 90%; (b) suboptimal input signal for q = 0.05 and Deff = 90%; (c) suboptimal input signal for q = 0.10 and Deff = 90%; and (d) suboptimal input signal for q = 0.40 and Deff = 90%.
Entropy 20 00528 g001
Figure 2. The system parameter identification block diagram.
Figure 2. The system parameter identification block diagram.
Entropy 20 00528 g002
Figure 3. Ellipsoidal confidence regions of LS parameter estimates (a); and LEL parameter estimates for different inputs: D-optimal input signal (black dotted, Deff = 90%, q = 1 × 10−6), Suboptimal input signal (green dashed, Deff = 90%, q = 0.10), Suboptimal input signal (red dash-dotted, Deff = 90%, q = 0.40), and Step input signal (blue solid line) (b).
Figure 3. Ellipsoidal confidence regions of LS parameter estimates (a); and LEL parameter estimates for different inputs: D-optimal input signal (black dotted, Deff = 90%, q = 1 × 10−6), Suboptimal input signal (green dashed, Deff = 90%, q = 0.10), Suboptimal input signal (red dash-dotted, Deff = 90%, q = 0.40), and Step input signal (blue solid line) (b).
Entropy 20 00528 g003
Figure 4. The comparison of the estimates values obtained using LEL (black solid line), and LS (red dashed line) estimators for various initial conditions, different noise variance, and D-efficiency value 90%: (a) confrontation between rival estimators for parameter estimate a based on the 80 runs; and (b) confrontation between rival estimators for parameter estimate b based on the 80 runs.
Figure 4. The comparison of the estimates values obtained using LEL (black solid line), and LS (red dashed line) estimators for various initial conditions, different noise variance, and D-efficiency value 90%: (a) confrontation between rival estimators for parameter estimate a based on the 80 runs; and (b) confrontation between rival estimators for parameter estimate b based on the 80 runs.
Entropy 20 00528 g004
Table 1. Comparison of the performance index components for the Deff = 100%.
Table 1. Comparison of the performance index components for the Deff = 100%.
Deff/DoptFIMJ1qqJ2tf [s]
100%−43.871.001 × 10−61 × 10−410.00
1.060.054.7910.57
1.100.108.9911.00
1.130.2017.4211.36
1.150.3025.9111.53
1.160.4034.4411.62
1.170.5042.9711.70
Table 2. Comparison of the performance index components for the Deff = 90%.
Table 2. Comparison of the performance index components for the Deff = 90%.
Deff/DoptFIMJ1qqJ2tf [s]
90%−35.600.8801 × 10−61 × 10−48.79
0.9230.055.009.23
0.9670.108.709.67
0.9910.2016.869.91
1.0410.3025.1410.04
1.0110.4033.4410.11
Table 3. Comparison of the performance index components for the Deff = 80%.
Table 3. Comparison of the performance index components for the Deff = 80%.
Deff/DoptFIMJ1qqJ2tf [s]
80%−28.100.7631 × 10−61 × 10−47.63
0.8070.054.358.07
0.8370.108.278.37
0.8580.2016.258.58
0.8660.3024.228.66
Table 4. The percentage accuracy rates for different values of q, and Deff = 90%.
Table 4. The percentage accuracy rates for different values of q, and Deff = 90%.
Indexq = 1 × 10−6q = 0.10
EstimatorLSLELLSLEL
Parametersabababab
average value [%]9.088.996.016.5710.9711.758.058.68
maximum value [%]61.2763.5447.5440.4455.6957.8240.7442.18
minimum value [%]2.0 × 10−21.6 × 10−12.1 × 10−43.3 × 10−41.9 × 10−17.0 × 10−21.5 × 10−41.2 × 10−4
Table 5. The percentage accuracy rates for different values of q, and Deff = 90%.
Table 5. The percentage accuracy rates for different values of q, and Deff = 90%.
Indexq = 0.40Step Input Signal
EstimatorLSLELLSLEL
Parametersabababab
average value [%]12.2912.459.029.4914.3514.239.8310.56
maximum value [%]63.2747.4349.9843.1583.1986.9554.4959.76
minimum value [%]4.6 × 10−29.6 × 10−11.3 × 10−69.5 × 10−51.9 × 10−12.3 × 10−24.0 × 10−41.5 × 10−4

Share and Cite

MDPI and ACS Style

Jakowluk, W. Free Final Time Input Design Problem for Robust Entropy-Like System Parameter Estimation. Entropy 2018, 20, 528. https://doi.org/10.3390/e20070528

AMA Style

Jakowluk W. Free Final Time Input Design Problem for Robust Entropy-Like System Parameter Estimation. Entropy. 2018; 20(7):528. https://doi.org/10.3390/e20070528

Chicago/Turabian Style

Jakowluk, Wiktor. 2018. "Free Final Time Input Design Problem for Robust Entropy-Like System Parameter Estimation" Entropy 20, no. 7: 528. https://doi.org/10.3390/e20070528

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop