Stochastic distribution tracking control for stochastic non-linear systems via probability density function vectorisation

Abstract

This paper presents a new control strategy for stochastic distribution shape tracking regarding non-Gaussian stochastic non-linear systems. The objective can be summarised as adjusting the probability density function (PDF) of the system output to any given desired distribution. In order to achieve this objective, the system output PDF has first been formulated analytically, which is time-variant. Then, the PDF vectorisation has been implemented to simplify the model description. Using the vector-based representation, the system identification and control design have been performed to achieve the PDF tracking. In practice, the PDF evolution is difficult to implement in real-time, thus a data-driven extension has also been discussed in this paper, where the vector-based model can be obtained using kernel density estimation (KDE) with the real-time data. Furthermore, the stability of the presented control design has been analysed, which is validated by a numerical example. As an extension, the multi-output stochastic systems have also been discussed for joint PDF tracking using the proposed algorithm, and the perspectives of advanced controller have been discussed. The main contribution of this paper is to propose: (1) a new sampling-based PDF transformation to reduce the modelling complexity, (2) a data-driven approach for online implementation without model pre-training, and (3) a feasible framework to integrate the existing control methods.

Introduction

The stochastic distribution control was presented for a class of non-Gaussian stochastic systems in the late 1990s (Wang, 1999), where the randomness of the system output can be controlled by adjusting the shape of the output probability density function (PDF). As an important research topic, stochastic distribution control inspires other topics such as the fault diagnosis in non-Gaussian systems (Guo and Wang, 2005; Yao et al., 2012), networked Direct Current (DC) motor control (Ren et al., 2015), probabilistic decoupling (Zhang et al., 2017), performance enhancement (Zhou et al., 2016), data-based identification (Zhang and Sepulveda, 2017), non-Gaussian filtering (Zhang and Yin, 2018; Zhao and Mili, 2017), operational control (Ding et al., 2012; Zhang and Hu, 2018), multi-path estimation (Cheng et al., 2018), industry 4.0 (Trovati et al., 2019), and so forth. In practice, tracking the given desired PDF is required in many process control and manufacturing processes, such as the quality control for paper-making (Wang, 1998).

A key step in PDF control is to establish the dynamic system model. Most of the existing results can be divided into two groups (Ren et al., 2019): (1) establish the PDF of the system output analytically if the stochastic distribution of the system noise is known (Guo and Wang, 2010). Following the analytical formulation, the full-probabilistic design can be achieved (Zhou and Herzallah, 2020) and the non-Gaussian filtering problem was investigated using the same approach (Guo and Wang, 2006). However, the PDF evolution is difficult to obtain considering the complex nature of stochastic systems and that the stochastic distribution of the system noise is normally unknown in practice. (2) represent the PDF as a weighing vector of the neural network, such as B-spline neural networks, where the dynamics of the system output can be reflected by the time-variant weighting vector (Wang, 2012). Although the PDF of the system noise is not necessarily known, the pre-training of the neural network can be applied to approximate the real-time PDF with the weighting-based representation. One problem of this approach is that it leads to a time delay for the control system design.

To overcome the aforementioned problems, it is important to develop a new approach combining the features of the two approaches above, where the new approach can be implemented in simple formulation and without requiring pre-training. This motivates the model development in this paper, where the evolutionary method is used to obtain the time-variant PDF. However, the analytically PDF will not be used for control design due to the computational cost. Alternatively, the sampling-based PDF vectorisation method is introduced, which converts the PDF to a vector with reduced dimension. Thus, the vector-based PDF model can be used in control system design. Note that the pre-training is not required for obtaining the vector, which makes the main difference from the neural network approach. Basically, the desired PDF can be converted into a vector. Then, the PDF tracking problem is re-written as a vector assignment problem. Any existing control methods can be used to minimise or eliminate the error between the reference vector and the time-variant vector along the time horizon. After identifying the model of the PDF dynamics, in this paper, we use Proportional-integral-derivative (PID) as a standard design, which is widely used in industry. Moreover, the stability of the proposed method is analysed to guarantee the tracking performance

In addition, the data-driven approach is also discussed as an extension of the presented strategy. The output PDF can be estimated by a sliding window that is known as the kernel density estimation (KDE) (Odiowei and Cao, 2009; Tang et al. 2020). Then, the data-based estimated PDF can be used to replace the evolutionary PDF. Note that the data-driven vectorisation of PDF can also be attained via Monte Carlo methods (Zhang and Wang, 2020), which can be considered as the PDF discretisation with zero-order hold. In addition, a first-order B-spline method (Wang, 2012) will lead to the equivalent histogram compared with Monte Carlo methods (Zhang and Wang, 2020). However, the modelling accuracy of these methods depends on the selection of the intervals in sample space, for which the pre-specified intervals may not sensitively reflects the dynamics of the PDFs. The hybrid modelling can be further discussed as the future perspective of the PDF tracking research where the model-data fusion is inspired for industrial process control and optimisation.

In practice, the data-driven approach will provide more flexibilities in terms of the implementation as the analytical formula of the stochastic distribution is difficult to obtain. In addition, the data-driven approach will be naturally included into the artificial intelligence design, where many recent design options are available. All these benefits can be considered in the future work.

The remainder of this paper is organised as follows. In the next section, preliminaries and the PDF representation are given. After that, the main result on PID control design is presented. Following the previous section, the stability analysis is given. The data-driven extension of the presented algorithm is proposed in the subsequent section. Then, numerical simulation results are demonstrated to validate the proposed method. In addition, the multi-output systems and advanced controller design are discussed at the end as an extension to the presented framework. Conclusions are given in the final part of the paper.

Preliminaries and PDF evolution

Consider the following general stochastic non-linear system with single input and single output

$$\begin{array}{c}{x}_{k+1}=f(x_k,u_k)+w_k\hfill \\ y_k=h\left(x_k\right)+v_k\hfill \end{array}$$

(1)

where $x\in {\mathbb{R}}^n$ , $y\in {\mathbb{R}}^1$ and $u\in {\mathbb{R}}^1$ stand for the system state vector, system output and control input, respectively. $w\in {\mathbb{R}}^n$ and $v\in {\mathbb{R}}^1$ are given as the random noises where the stochastic distributions are known. $k$ is the sampling index. $f:{\mathbb{R}}^n\times {\mathbb{R}}^1\to {\mathbb{R}}^n$ and $h:{\mathbb{R}}^n\to {\mathbb{R}}^1$ are non-linear functions. The following assumptions are specified for the investigated system (1):

Assumption 1: There exist two positive real numbers $L_1$ and $L_2$ , such that, the following inequality always holds for $k=1,2,\dots$ .

$$||f(x_k,u_k)-f(x_{k-1},u_{k-1})||⩽L_1||\mathrm{\Delta }{x}_k||+L_2||\mathrm{\Delta }{u}_k||$$

where $\mathrm{\Delta }{x}_k=x_k-x_{k-1}$ and $\mathrm{\Delta }{u}_k=u_k-u_{k-1}$ .

Assumption 2: The non-linear function $h(·)$ meets Lipschitz condition where exists a real positive number $L_3$ , such that

$$||h\left(x_k\right)-h\left(x_{k-1}\right)||⩽L_3||\mathrm{\Delta }{x}_k||$$

Assumption 3: The Jacobian determinant ${\mathrm{\Xi }}_k$ satisfies

$${\mathrm{\Xi }}_k=\left|\det \frac{\partial v}{\partial {\tau }_2}\right|\ne 0$$

where $v=H^{-1}(x_{k-1},u_{k-1},{\tau }_1,{\tau }_2)$ and $H(·)$ is a non-linear function. In particular

$${\tau }_2=H(x_{k-1},u_{k-1},{\tau }_1,v)$$

while ${\tau }_1$ , ${\tau }_2$ , $v$ denote the random variables of $w_k$ , $y_k$ and $v_k$ .

To describe the evolution of the output PDF for stochastic non-linear system (1), the following lemma has been recalled (Yin and Guo, 2012; Zhou et al., 2017).

Lemma 1: For the stochastic system (1) with Assumption 3, the PDF of the system output is formulated by

$$\begin{array}{c}{\gamma }_k\left({\tau }_2\right)={\int }_{\mathrm{\Omega }}{\gamma }_v\left(H^{-1}(x_{k-1},u_{k-1},{\tau }_1,{\tau }_2)\right){\gamma }_w\left({\tau }_1\right)\hfill \\ \times {\left|\det \frac{\partial H^{-1}(x_{k-1},u_{k-1},{\tau }_1,{\tau }_2)}{\partial {\tau }_2}\right|}^{-1}d{\tau }_1\hfill \end{array}$$

(2)

where ${\gamma }_v\left(v\right)$ and ${\gamma }_w\left({\tau }_1\right)$ stand for the PDFs of $v_k$ and $w_k$ . $\mathrm{\Omega }$ denotes the sample space of the system output.

Using the lemma above, the PDF of system output ${\gamma }_k({\tau }_1,u_k)$ can be obtained for each time instant $k$ . The PDF problem can be described as follows

$$min\underset{k\to \infty }{lim}{\int }_{\mathrm{\Omega }}{({\gamma }_{\mathrm{ref}}-{\gamma }_k({\tau }_1,u_k))}^{2}d{\tau }_1$$

(3)

where ${\gamma }_{\mathrm{ref}}$ denotes the desired PDF of the system output. Equation (3) will approach zero asymptotically if the perfect tracking is achieved. To simplify the computation in (3), we can define a set of base points in the sample space for output, such as ${\alpha }_1,{\alpha }_2,\dots ,{\alpha }_m$ where $m⩾n$ is a positive integer. Then, at time k, the PDF can be vectorised using the sampling approach as follows

$$z_k={\left[\begin{array}{cccc}{\gamma }_k\left({\alpha }_1\right)\hfill & {\gamma }_k\left({\alpha }_2\right)\hfill & \dots \hfill & {\gamma }_k\left({\alpha }_m\right)\hfill \end{array}\right]}^T$$

(4)

where $z_k\in {\mathbb{R}}^m$ is formulated as a vector that is defined as the vector-valued PDF representative state. Notice that ${\gamma }_k$ can be approximated by $z_k$ when the dimension of $z_k$ is large enough. To demonstrate the relationship of the base points and the PDF re-establishment, Figure 1 illustrates the samplings of the PDF.

Figure 1.

PDF representation using vectorisation which can be considered as PDF sampling.

Remark 1: The selection of the vector points in (4), including the number and locations, will influence the accuracy of PDF approximation. In principle, more points will lead to closer approximation but with higher computational load. In practice, we hope to balance the modelling accuracy and computational cost. Among the $m$ points in the vector, we include all local maximum and minimum points in output PDF.

The desired PDF ${\gamma }_{\mathrm{ref}}$ can also be converted into a vector with dimension $m$ , which is written as $z_{\mathrm{ref}}$ . As a result, control objective (3) can be re-written using the vector-valued PDF representative state to replace the original continues PDF, which can be formulated as follows

$$\underset{k\to \infty }{lim}{\tilde{z}}_k\to 0$$

(5)

where ${\tilde{z}}_k=z_{\mathrm{ref}}-z_k$ is defined as the error of vector assignment at time k.

In particular, the selection of minimum set of the base points can be implemented using the inflection points of the desired PDF. For each two inflection points of the PDF, the curve between them is positive and monotonous. As the integral of the PDF over the definition space is equal to one, the inflection points can uniquely represent the shape of the desired PDF. Thus, the minimum dimension of the vector-valued PDF representative state can be obtained. However, the minimum set of vector points cannot reflect the transient part of the PDF dynamics. Additional points can be pre-selected between each two inflection points, which would increase the accuracy of the control performance when additional transient information has been considered into the PDF control.

Remark 2: Notice that $z_k$ and $z_{\mathrm{ref}}$ can be obtained directly from the PDF formula using Lemma 1. This means that the vector can be used in real time without pre-training.

PID controller design algorithm

Using the PDF vectorisation, a new output PDF model can be established where many existing controller design approaches can be applied such as PID, Linear-quadratic regulator (LQR), and so forth. For instance, a PID controller with the following structure

$$u_k=K_P{\tilde{z}}_k+K_I\sum _{i=1}^k{\tilde{z}}_i+K_D({\tilde{z}}_k-{\tilde{z}}_{k-1})$$

(6)

where $K_P\in {\mathbb{R}}^m$ , $K_I\in {\mathbb{R}}^m$ and $K_D\in {\mathbb{R}}^m$ stand for the proportional gain, integral gain and derivative gain, respectively.

The parameters of the controller can be determined by various methods; in this paper, we apply the time-domain design approach using linear matrix inequality (LMI). Firstly, the model should be built up to describe the dynamics of the PDF based upon the transformed vector-valued PDF representative state. In particular, the following linear format is used

$$z_{k+1}=A{z}_k+B{u}_k$$

(7)

where $A\in {\mathbb{R}}^{m\times m}$ and $B\in {\mathbb{R}}^{m\times 1}$ stand for the coefficient matrices of the system model. Then, the model (7) can be further represented as

$$z_{k+1}=\mathrm{\Theta }{\left[\begin{array}{cc}{z}_k^T\hfill & u_k^T\hfill \end{array}\right]}^T$$

(8)

where the entries in $\mathrm{\Theta }=\left[\begin{array}{cc}A\hfill & B\hfill \end{array}\right]\in {\mathbb{R}}^{m\times (m+1)}$ can be identified using the least square algorithm. Note that $z_i$ , $u_i$ with $\forall i<k$ are available for each time instant $k$ .

As a result, we can further rewrite the PID formula as follows

$$\begin{array}{c}{u}_k=K{\left[\begin{array}{cccc}{\tilde{z}}_k^T\hfill & p_k^T\hfill & {\tilde{z}}_k^T\hfill & -{\tilde{z}}_{k-1}^T\hfill \end{array}\right]}^T\hfill \\ p_k=p_{k-1}+{\tilde{z}}_k\hfill \end{array}$$

(9)

where $p_k=\sum _{i=1}^k{\tilde{z}}_i$ and $K=\left[\begin{array}{ccc}{K}_P\hfill & K_I\hfill & K_D\hfill \end{array}\right]$ .

Defining $\overline{z}={\left[\begin{array}{ccc}{z}_k^T\hfill & p_k^T\hfill & z_{k-1}^T\hfill \end{array}\right]}^T$ as the generalised vector of the system output vector, the closed-loop system model can be given by

$$\begin{array}{c}{\overline{z}}_{k+1}=(\overline{A}+\overline{B}K\overline{C}){\overline{z}}_k+{\left[\begin{array}{ccc}B{K}_P\hfill & I\hfill & 0\hfill \end{array}\right]}^T{z}_{\mathrm{ref}}\hfill \\ z_k=\left[\begin{array}{ccc}I\hfill & 0\hfill & 0\hfill \end{array}\right]{\overline{z}}_k\hfill \end{array}$$

(10)

where $I$ stands for identity matrix and

$$\overline{A}=\left[\begin{array}{ccc}A\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-I}& I\hfill & 0\hfill \\ \multicolumn{1}{c}{I}& 0\hfill & 0\hfill \end{array}\right],\overline{B}=\left[\begin{array}{c}B\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right],\overline{C}=\left[\begin{array}{ccc}-I\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& I\hfill & 0\hfill \\ \multicolumn{1}{c}{-I}& 0\hfill & I\hfill \end{array}\right]$$

(11)

Based on the identified linear model, the parameter matrix $K$ can be selected and the result has been summarised in the following theorem.

Theorem 1: There exists a gain matrix $K$ that makes the dynamics of the vector-based linear model (7) asymptotically stable with the extended coefficient matrices (10) and (11), such that the parametric matrix $K=W^{-1}Y$ satisfies the following LMI

$$\left[\begin{array}{cc}-M\hfill & M\overline{A}+\overline{B}Y\overline{C}\hfill \\ \multicolumn{1}{c}{{(M\overline{A}+\overline{B}Y\overline{C})}^T}& (\beta -1)M\hfill \end{array}\right]<0$$

(12)

where $M$ stands for a symmetric positive definite matrix, $M\overline{B}=\overline{B}W$ and $\beta ⩾0$ denotes the decay rate.

Proof: The following Lyapunov function candidate can be selected for ensuring the stability of the control design

$$V_k={\overline{z}}_k^{T}M{\overline{z}}_k$$

(13)

where $M$ stands for the real positive definite symmetric matrix $M>0$ .

Based upon the presented Lyapunov function candidate, we have the following equation

$$V_{k+1}-V_k={\overline{z}}_k^T({(\overline{A}+\overline{B}K\overline{C})}^{T}M(\overline{A}+\overline{B}K\overline{C})-M){\overline{z}}_k$$

(14)

We can further consider a decay rate for the Lyapunov function candidate such that the following inequality holds

$$V_{k+1}-V_k⩽-\beta V_k$$

(15)

where $0⩽\beta <1$ .

Thus, the condition of system stability can be described as follows

$${(\overline{A}+\overline{B}K\overline{C})}^{T}M(\overline{A}+\overline{B}K\overline{C})-(1-\beta )M<0$$

(16)

To rewritten the inequality into LMI, we can further introduce $Y=\mathrm{WK}$ and $M\overline{B}=\overline{B}W$ , thus the LMI (12) is obtained using Schur complement and the PID gain can be selected as $K=W^{-1}Y$ , which ends the proof.

Remark 3: With the vector-based modelling, the computation has been simplified, which potentially benefits implementation to applications.

Stability

In this section, the stability of the proposed output PDF control system is analysed. Combining the system model (1) and PID law (6), the following equations are obtained

$$\mathrm{\Delta }{x}_{k+1}=f(x_k,u_k)-f(x_{k-1},u_{k-1})+w_k-w_{k-1}$$

(17)

and

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Delta }{u}_k=(K_P+K_D)({\tilde{z}}_k-{\tilde{z}}_{k-1})}\\ -K_D({\tilde{z}}_{k-1}-{\tilde{z}}_{k-2})+K_I{\tilde{z}}_k\hfill \end{array}$$

(18)

where $\mathrm{\Delta }{x}_{k+1}=x_{k+1}-x_k$ and $\mathrm{\Delta }{u}_k=u_k-u_{k-1}$ .

Based upon the two assumptions of $f(·)$ , equation (17) results in

$$||\mathrm{\Delta }{x}_{k+1}||⩽L_1||\mathrm{\Delta }{x}_k||+L_2||\mathrm{\Delta }{u}_k||+||w_k-w_{k-1}||$$

(19)

Notice that there always exist two real numbers $\underset{¯}{\lambda }$ and $\overline{\lambda }$ , such that the following inequality holds

$$\underset{¯}{\lambda }{\tilde{z}}_k⩽{\tilde{z}}_{k-j}-{\tilde{z}}_{k-j-1}⩽\overline{\lambda }{\tilde{z}}_k,j=0,1,\dots ,k-1$$

(20)

where $\underset{¯}{\lambda }$ and $\overline{\lambda }$ stand for the upper bound and lower bound coefficients.

Equation (18) leads to the following result using the norm operation

$$||\mathrm{\Delta }{u}_k||⩽(\overline{\lambda }(K_P+K_D)-\underset{¯}{\lambda }{K}_D+K_I)||{\tilde{z}}_k||$$

(21)

Since the vector $z_k$ will track the reference under the control input, there always exists a positive real number $\overline{\sigma }>0$ , such that the following inequality holds

$$||{\tilde{z}}_k||⩽\overline{\sigma }||{\tilde{x}}_k||$$

(22)

Calculating the mean-value of equation (19) and substituting equations (21)–(22) to equation (19), we can have

$$E\{||\mathrm{\Delta }{x}_{k+1}||\}⩽\mathrm{\Xi }E\{||\mathrm{\Delta }{x}_k||\}+{\delta }_w$$

(23)

where

$$\mathrm{\Xi }=||L_1+\overline{\sigma }{L}_2(\overline{\lambda }(K_P+K_D)-\underset{¯}{\lambda }{K}_D+K_I)||⩽1$$

(24)

as ${\delta }_w$ denotes the upper bound of $E\{||w_k-w_{k-1}||\}$ .

Using the assumption of function $h(·)$ , the system output $y$ is bounded at all time.

As a summary of the analysis in this section, a theorem is given to indicate the stability condition of the proposed output PDF control algorithm.

Theorem 2: The system output of the stochastic non-linear system (1) is bounded in mean-norm sense using the PID control design (6) with the system Assumptions 1 and 2, if the parameters of the controller can be selected to make the positive coefficient $\mathrm{\Xi }⩽1$ .

Proof: The proof of this theorem has been illustrated above.

Data-driven extension with KDE

As aforementioned in the Introduction, the KDE can be adopted to estimate the PDF using the data of system output. In this case, the evolutionary PDF can be replaced by the estimated PDF that further reduces the complexity of the computing. Since the collected system output data can be denoted as $y_1,y_2,\dots ,y_k$ , the PDF can be approximated as follows

$${\hat{\gamma }}_k\left({\tau }_1\right)=\frac{1}{k}\sum _{i=1}^{k}G({\tau }_1-y_i)$$

(25)

where $G(·)$ stands for the Gaussian kernel and ${\hat{\gamma }}_k\left({\tau }_1\right)$ denotes the output PDF for each time instant $k$ . As a result, the vector $z_k$ can be re-written as an estimated vector, that is

$${\hat{z}}_k={\left[\begin{array}{cccc}{\hat{\gamma }}_k\left({\alpha }_1\right)\hfill & {\hat{\gamma }}_k\left({\alpha }_2\right)\hfill & \dots \hfill & {\hat{\gamma }}_k\left({\alpha }_m\right)\hfill \end{array}\right]}^T$$

(26)

Note that the estimated vector ${\hat{z}}_k$ is a close approximation to the true value of the vector from PDF vectorisation if the data is sufficient. The estimation will not affect the stability analysis above. The estimation error can be merged into the identification error, which will be handled by the robustness of the controller design.

Based on the stability analysis, the flowchart has been produced in Figure 2 to demonstrate the procedure of the presented algorithm, where the flow with hollow arrows indicates the data-based extension of the presented framework. In particular, the framework has been shown to include three components: PDF evolution, PDF dynamics modelling and PDF tracking control.

Figure 2.

The flowchart of the presented PDF tracking control framework using PID controller.

Simulation

To illustrate the effectiveness of the proposed stochastic distribution control algorithm, a numerical example is demonstrated in this section and the system model is formulated as follows

$$\begin{array}{c}{x}_{k+1}=x_k\cos (x_k+0.03)+0.8u_k+w_k\hfill \\ \multicolumn{1}{c}{y_k=0.1x_k+v_k}\end{array}$$

where $w_k$ and $v_k$ stand for the process noise and measurement noise, respectively. Both $w_k$ and $v_k$ are zero-mean Gaussian noises and the variances of them are equal to $1$ and $0.1$ . Since the non-linear term exists in the model, the system output can be considered as a non-Gaussian random variable. The non-linear dynamics will re-shape the PDF of the system output even if the system noises are subjected to Gaussian distribution.

Suppose that the desired PDF is given as Gamma distribution ${\gamma }_{\mathrm{ref}}(x,a,b)=\frac{1}{b^a(a-1)!}{x}^{(a-1)}{e}^{\frac{-x}{b}}$ , where $a=1$ and $b=2$ . Therefore, the reference vector-valued PDF representative state can be obtained as $z_{\mathrm{ref}}=[0.005,0.2353,0.1570,0.0027]$ where the base pointes are pre-specified as $[-5,-1.667,1.667,5]$ in the output sample space. The parameter $K$ of the PID controller is selected as $K_P=\mathrm{diag}\{2000,2000,2000,2000\}$ , $K_I=\mathrm{diag}\{250,250,250,250\}$ and $K_D=0$ . The following figures are obtained to show the simulation results. The PID controller will take action for the closed-loop system since $k⩾100$ . Before that, the system output should be simulated for data collection, while the KDE is used to estimate the PDF. The collected data can also be used for vector-based dynamic model (7) identification if the controller is model-based design; for example, LQR, and so forth.

Figure 3 shows the vectorised PDF $z_k$ at $k=100$ and $k=400$ . It has been seen that the vector $z_k$ converges to $z_{\mathrm{ref}}$ . The control input with PID controller is shown in Figure 4, where the control input signal is bounded and convergent. In Figure 5, the 3D mesh of vector has been given to show the dynamics of the vector-valued PDF representative state where the vector $z_k$ is changed by the control input. Based upon $z_k$ and the interpolation using Matlab, the pseudo PDFs of the system output $y$ are shown in Figure 6. In order to re-store the real-time PDF, the interpolation approach has been adopted here. To illustrate the evolution of the PDF with details, Figure 7 is given where the PDFs at various time instants are shown. Note that the PDF converges to the desired target PDF along the time horizon. Furthermore, Figure 8 shows the measured system output along $k$ . It has been shown that the system output is also bounded with PID design, which has validated the result of the stability analysis. Moreover, the tracking error vector of the PDF representative state $\tilde{z}$ is indicated by Figure 9. Considering the mean-square criterion of the tracking error $\tilde{z}$ , the curves in Figure 10 show the attenuation of the tracking error. It also implies the tracking error convergence of the presented PDF control algorithm. In addition, compared with Figure 6, one 3D mesh of the evolutionary PDF of the system output is shown in Figure 11, where the control input signal $u_k$ has been fixed as $u_k=0$ for all $k$ and the system performance illustrates that the system output PDF of the investigated stochastic system is not able to achieve the desired tracking performance without control input.

Figure 3.

The vector-valued PDF representative state $z$ of the system output $y$ comparing with $k=100,1000$ and $z_{\mathrm{ref}}$ .

Figure 4.

The PID control input for the investigated closed-loop system.

Figure 5.

The 3D mesh for the vector $z$ of the system output $y$ while $z$ converges to $z_{\mathrm{ref}}$ along $k$ .

Figure 6.

The pseudo PDF of the system output $y$ based on the curve-fitting, PDF representative vector and data interpolation.

Figure 7.

The interpolation-based pseudo PDF of the system output $y$ at various time instants.

Figure 8.

The system output $y$ of the stochastic system while $y$ is bounded based on the presented data-driven design.

Figure 9.

The error of the probability density states $\tilde{z}$ using presented data-driven approach.

Figure 10.

The mean-square error of the probability density states $\tilde{z}$ using presented data-driven approach.

Figure 11.

The vector-valued PDF representative state of the system output $y$ without control input where $u=0$ which shows that the shape of the system output PDF will not be changed along $k$ .

It has been shown from the numerical study that the presented algorithm is convenient for implementing without pre-training of weights. Using the data-driven approach, the PDF evolution can be replaced by a vector dynamics while the system identification may not be necessary as the PID parameters turning is achievable by trial and error. We can consider the presented framework as a model-free design for PDF tracking problem.

Further discussion

Multi-output systems extension

In this paper, we only investigate the single-input, single-output (SISO) systems where the PDF of the system output can be visualised by the 2D curve at each sampling instant. For complex industrial processes, the multi-input, multi-output (MIMO) systems should be discussed. The multi-variable system output will lead to a multi-dimensional PDF which is a joint PDF. Thus, the PDF tracking problem for multi-output systems can be summarised as a joint PDF tracking problem.

Suppose that the system outputs are $\overline{n}$ -dimensional vector, $y_i\in {\mathbb{R}}^{\overline{n}},i=1,2,\cdots ,\overline{n}$ , subject to random noises. As a result, the joint PDF can be denoted as $\gamma ({\omega }_1,{\omega }_2,\cdots ,{\omega }_{\overline{n}})$ where ${\omega }_i,i=1,2,\cdots ,\overline{n}$ stand for the random variable of each system output.

The sets of base points for each system output will be pre-specified as follows

$${\mathrm{\Omega }}_{y_i}=[{\alpha }_{i,1},{\alpha }_{i,2},\dots ,{\alpha }_{i,m_i}],i=1,2,\cdots$$

(27)

where $m_i$ denotes the dimension of vector-valued PDF representative state for each system output. For each output variable, $m_i$ can be selected individually as a large positive integer. It means that the vector size can be different for each output variable. In practice, we can simply assume that $m_1=m_2=\cdots =m_{\overline{n}}$ to reduce the complexity of the expression without loss of the generality.

The multi-dimensional vector-valued PDF representative state can be obtained by substituting the base points into the joint PDF, then the joint PDF of the system output can be represented by a single vector. Take an example with two system outputs $y_1$ and $y_2$ , suppose that $m_1=m_2$ , thus we have the following square matrix to represent the joint PDF $\gamma ({\omega }_1,{\omega }_2)$

$$\left[\begin{array}{cccc}\gamma ({\alpha }_{1,1},{\alpha }_{2,1})\hfill & \gamma ({\alpha }_{1,1},{\alpha }_{2,2})\hfill & \cdots \hfill & \gamma ({\alpha }_{1,1},{\alpha }_{2,m_2})\hfill \\ \multicolumn{1}{c}{\gamma ({\alpha }_{1,2},{\alpha }_{2,1})}& \gamma ({\alpha }_{1,2},{\alpha }_{2,2})\hfill & \cdots \hfill & \gamma ({\alpha }_{1,2},{\alpha }_{2,m_2})\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\gamma ({\alpha }_{1,m_1},{\alpha }_{2,1})}& \gamma ({\alpha }_{1,m_1},{\alpha }_{2,2})\hfill & \cdots \hfill & \gamma ({\alpha }_{1,m_1},{\alpha }_{2,m_2})\hfill \end{array}\right]$$

(28)

Vectorising the matrix, the vector-valued PDF representative state for the joint PDF can be formulated as a vector. In particular, we can choose each row of the matrix and connect them as a vector, alternatively, the column approach can be used to form the vector. Both types of vectorisation will give the same dynamics of the joint PDF. For instance, the row-based vectorisation for the matrix above can be converted into the following vector

$$\begin{array}{c}z=[\gamma ({\alpha }_{1,1},{\alpha }_{2,1}),\gamma ({\alpha }_{1,1},{\alpha }_{2,2}),\cdots ,\gamma ({\alpha }_{1,1},{\alpha }_{2,m_2}),\hfill \\ \gamma ({\alpha }_{1,2},{\alpha }_{2,1}),\cdots ,\gamma ({\alpha }_{1,2},{\alpha }_{2,m_2}),\cdots ,\hfill \\ \gamma ({\alpha }_{1,m_1},{\alpha }_{2,1}),\cdots ,\gamma ({\alpha }_{1,m_1},{\alpha }_{2,m_2})]\hfill \end{array}$$

(29)

The multi-output system PDF tracking problem can be described uniformly using the proposed framework. In general, the vector can be denoted as $z\in {\mathbb{R}}^{\overline{m}},\overline{m}=\underset{i=1}{\overset{\overline{n}}{\Pi }}{m}_i$ while the coefficient matrices $A$ and $B$ are of the proper dimensions. Compared with the SISO systems, the MIMO systems will lead to a high-dimensional dynamic model. Then the parametric identification would become a challenge. Technically, the identified coefficient matrices would drop into ill-condition, where the additional operation is needed to analyse the condition number during the identification. In addition, the high-dimensional dynamic model will also lead to a complicated structure of the controller. The systematic synthesisation would be another challenge.

To implement the presented algorithm via a data-driven framework, the multi-dimensional KDE can be adopted, thus the joint PDF can be approximated using the collected data sets of the system outputs. The data-based estimated vector ${\hat{z}}_k$ for each sampling instant $k$ can be denoted using $\hat{\gamma }({\omega }_1,{\omega }_2,\dots ,{\omega }_{\overline{n}})$ . The data-driven approach is achievable once ${\hat{z}}_k$ is used for model identification and controller design.

Advanced controller design

Another extension of the presented framework can be investigated from the controller design point of view. Currently, we used PID as a standard controller for a linear system. However, the linear model may not represent complex output PDF systems, thus, non-linear models and un-modelled dynamics should be considered. Advanced controller design need to be developed on complex system models. In particular, the small-gain technique can be used to deal with the uncertainty and the robustness requirement (Ma et al., 2020). Fuzzy tracking control design also focuses on the uncertainties and robustness, where non-linear networked system model (Li and Park, 2018), high-order non-linear model (Zhao et al., 2015), event-triggered non-linear model (Li et al., 2019), non-strict feedback non-linear model (Wang et al., 2021) have been presented. All these mentioned methods can be implemented with the proposed data-driven approach.

Conclusions

This paper proposes an output PDF tracking control algorithm for a class of stochastic systems. Different from the existing neural network-based methods and PDF evolution solutions, the presented control algorithm converts the system output PDF to a vector-formed probability density states using vectorisation. The dynamics of the system output PDF is approximated by the time-variant vector-valued PDF representative state. Thus, the PDF tracking problem has been transformed into a vector assignment problem. To describe the relationship between the control input and the PDF representative vector, a vector-based model is established with linear dynamics assumed, where the vector is measurable using the PDF evolution and the PDF vectorisation. To achieve the PDF tracking, a LMI-based PID design is adopted to eliminate the distance between the reference vector and the real-time vector asymptotically. Following the theoretical analysis, the closed-loop stability of the system can be achieved where the PDF tracking error is convergent and the system output is bounded in the mean-norm sense.

To deal with the system noise with unknown distribution, the data-driven approach is further discussed where the vector can be estimated by KDE and the stability is still guaranteed. The numerical simulation results demonstrate the effectiveness of the presented design.

As the main contribution of the paper, a new PDF tracking control algorithm is proposed, which can be used in real-time control without pre-training of weights and potentially it can be implemented as a pure data-driven or even model-free algorithm. The computing complexity has been reduced by the vectorisation of continuous PDF, and existing controllers can be integrated into the framework.

Using the presented framework, the following three aspects can be considered as our future works: (1) data-driven optimisation using machine learning algorithms, (2) advanced controller design considering the un-modelled dynamics of PDF and (3) industrial applications with analysis and randomness attenuation.