A gaussian process regression method and apparatus for system identification

Abstract

The application provides a Gaussian process regression method and device for system identification, and belongs to the technical field of automatic control. The method comprises the following steps: obtaining an original data set of a model to be identified in an automatic control system; assigning an importance factor to each group of original data; updating the importance factor based on the Euclidean distance between input data; removing original data with low importance in the original data set to obtain a uniformly distributed data subset; grouping input data of the data subset into a training input matrix and grouping output data into a training output matrix; performing Gaussian process regression training on the data subset to obtain a covariance matrix considering data noise; and calculating an output prediction result corresponding to the input data to be predicted, so as to realize system identification. The application can extract a small sample set with uniform distribution, solve the problem of large calculation amount of Gaussian process regression when applied to system identification, and promote the realization of a high-precision stable controller.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of automatic control, and relates to a Gaussian process regression method and device for system identification. BACKGROUND

[0002] System identification is an important part of automatic control technology, and designing an adaptive control algorithm based on system identification can effectively improve the tracking performance and robustness of a control system. As a non-parametric regression method, Gaussian process regression (GPR) has attracted much attention from many scholars in the application of system identification. The basic principle of Gaussian process regression is to use Gaussian process prior and input-output data to obtain the output posterior distribution under the given input condition based on Bayesian inference. Compared with linear regression and neural network methods, Gaussian process regression has the advantages of adaptive hyperparameter acquisition, combination of prior knowledge, statistically meaningful output and strong generalization ability, and has advantages in high-dimensional, small-sample and nonlinear problems.

[0003] However, the calculation of Gaussian process regression requires matrix inversion, and there is a problem of large calculation amount when the training data is large. In order to reduce the calculation amount, the existing solutions mainly include data subset method, rank reduction approximation method and sparse pseudo-input method. Among them, the principle of the data subset method is to obtain a subset of the original data based on certain rules, and Gaussian process regression is only fitted to the data subset, thereby reducing the calculation amount. This method has better effect than the latter two in many cases. At present, the commonly used data subset method includes random selection and greedy algorithm. However, both of these two methods have limitations. Random selection can retain the distribution characteristics of the original data, but the data in the automatic control field is generally unevenly distributed, so the data subset obtained by random selection is also uneven, which reduces the fitting effect of the model. The application scope of the greedy algorithm is relatively small, and it often needs to be analyzed according to specific problems, so it cannot be effectively promoted in the field of automatic control. SUMMARY

[0004] The purpose of the present application is to overcome the shortcomings of the prior art and provide a Gaussian process regression method and device for system identification. The present application screens the data subset through the Euclidean distance between data, reduces the sample size while ensuring the effect of Gaussian process regression, thereby reducing the calculation amount, and has good adaptive ability for data sets of different dimensions and different distributions. For the unevenly distributed data commonly used in the field of automatic control, the present application can extract a small sample set with uniform distribution, thereby solving the problem of large calculation amount of Gaussian process regression in the application of system identification, and promoting the realization of high-precision stable controllers.

[0005] The embodiment of the first aspect of the present application provides a Gaussian process regression method for system identification, comprising:

[0006] An original data set of a model to be identified in an automatic control system is obtained, an importance factor is given to each group of original data, the importance factor is updated based on the Euclidean distance between input data in the original data, original data with low importance in the original data set is removed, and a uniformly distributed data subset is obtained;

[0007] Input data of the data subset is used to form a training input matrix, and output data of the data subset is used to form a training output matrix, Gaussian process regression training is performed on the data subset, and a covariance matrix considering data noise is obtained;

[0008] Based on the training output matrix, the training output matrix and the covariance matrix, an output prediction result corresponding to predicted input data in the model to be identified is calculated, so as to realize system identification.

[0009] In one specific embodiment of the present application, the uniformly distributed data subset is obtained, comprising:

[0010] 1) An original data set of a model to be identified in an automatic control system is obtained;

[0011] Wherein, the input dimension of the original data of the model to be identified is n i , and the output dimension is n o .

[0012] Let the original data set be S={(x (j) ,y (j) ),j=1,2,...,P}, wherein and are input data and output data of the jth group of data, respectively, and P is the total number of data groups;

[0013] 2) A data subset D is created, which is initially empty, and the capacity of the data subset D is C;

[0014] 3) Let j=1;

[0015] 4) Determine the number N of stored data groups in the data subset D:

[0016] If N=0, it indicates that the current data subset D is empty, the jth group of data (x (j) ,y (j) ) is added to D, the importance factor μ j of the group of data is set to 1, N is set to N+1, and the maximum value M i of the i th element in the first group to the jth group of input data in the original data set is recorded, wherein when j=1, i = 1, 2, …, n i ; set the minimum value of the i-th element in the 1st to j-th input data of the original data set as m i , where when j = 1, i = 1, 2, …, n i ; i = 1, 2, …, n i , k = 1, 2, …, j; add the current x (j) as the last added input data x last in the current data subset D, that is, x last = x (j) , and then go to step 6);

[0017] If N > 0, it means that the current data subset D is not an empty set, and the maximum and minimum values of each element in the 1st to j-th input data of the original data set are updated as follows:

[0018]

[0019] Calculate the Euclidean distance d between the input data x (j) in (x (j) , y (j) ) and the last added input data x last in the current data subset D:

[0020] d = ||x (j) -x last ||2 (2)

[0021] Let the set composed of the 1st to j-th input data of the original data set be X = {x (k) | k = 1, 2, …, j}, and calculate the volume V of the smallest n ni - lattice containing X in R i space:

[0022]

[0023] Calculate the receivable distance d a :

[0024]

[0025] Then go to step 5);

[0026] 5) Determine the Euclidean distance:

[0027] 5-1) If d ≤ d a , the j-th original data is not added to the current data subset D;

[0028] Then, determine j: if j = P, then all data in the original dataset has been traversed, and the current data subset D is the final data subset D; if j < P, then let j = j + 1, and then return to step 4);

[0029] 5-2) If d > d a Then the j-th group of original data (x) (j) ,y (j) Add it to the current subset D of data and set the importance factor μ for this set of data. j =1, let N = N+1, and add the currently added x (j) x is the last input data added to subset D. last That is, let x last =x (j) ;

[0030] Then update the Gaussian kernel parameter γ using the following formula when V > 0:

[0031]

[0032] Where s is the standard deviation of the Gaussian kernel; when V = 0, γ = 0;

[0033] Calculate the value of each data set (x) in the current subset D. (k) ,y (k) The importance of reducing gain G k k≠j:

[0034]

[0035] Update each data set (x) in the current subset D. (k) ,y (k) Importance factor μ k k≠j:

[0036] μ k :μ k ·G k (7)

[0037] Then proceed to step 6);

[0038] 6) Determine j:

[0039] If j = P, then all data in the original dataset has been traversed, and the current data subset D is the final data subset D; if j < P, then N is judged: if N > C, it means that the data subset is full, and the data with the smallest importance factor (x) in the current data subset D is selected. (r) ,y (r) ), Move out D, then let j = j + 1, and then return to step 4); if N ≤ C, let j = j + 1, and then return to step 4).

[0040] In one specific embodiment of the present application, the covariance matrix considering data noise comprises:

[0041] Taking the data subset D as the training set, the number of data groups of the training set is N, the training input matrix of the training set is X train = [x (1) , x (2) ,..., x (N) ] T , and the training output matrix is Y train = [y (1) , y (2) ,..., y (N) ] T , and are column vectors, respectively representing the i-th input data and the i-th output data in the data subset D, i = 1, 2,..., N;

[0042] The square exponential kernel is used to calculate the covariance matrix K y :

[0043]

[0044] wherein K y is the covariance matrix considering data noise, K is the covariance matrix without considering data noise, is the noise variance of the output data, I is an N × N unit matrix, and p is the square exponential kernel parameter.

[0045] In one specific embodiment of the present application, the calculation of the output prediction result corresponding to the input data to be predicted based on the training output matrix, the training output matrix and the covariance matrix comprises:

[0046] Let the input data to be predicted be j = 1, 2,..., n, and n is the number of input data to be predicted;

[0047] The covariance matrix K * between the training data and the prediction data and the covariance matrix K ** between the prediction data are calculated:

[0048]

[0049] Then the output prediction matrix Y * and the covariance Σ * are calculated:

[0050]

[0051] wherein, is a prediction value corresponding to a covariance matrix Σ * the jth element on the diagonal is a variance estimate of the prediction value

[0052] The second aspect of the present application provides a Gaussian process regression device for system identification, comprising:

[0053] A data subset construction module is configured to obtain original data groups of a model to be identified in an automatic control system to form an original data set, assign a corresponding importance factor to each original data group, update the importance factor based on the Euclidean distance between input data in the original data, eliminate original data with low importance in the original data set, and obtain a uniformly distributed data subset.

[0054] A Gaussian process regression training module is configured to group input data of the data subset into a training input matrix, group output data of the data subset into a training output matrix, and obtain a covariance matrix considering data noise by performing Gaussian process regression training on the data subset.

[0055] A prediction module is configured to calculate an output prediction result corresponding to input data to be predicted in the model to be identified based on the training output matrix, the training output matrix, and the covariance matrix, so as to realize system identification.

[0056] The third aspect of the present application provides an electronic device, comprising:

[0057] at least one processor; and a memory connected to the at least one processor in communication;

[0058] wherein the memory stores instructions executable by the at least one processor, and the instructions are configured to perform the above-mentioned Gaussian process regression method for system identification.

[0059] The fourth aspect of the present application provides a computer readable storage medium, which stores computer instructions for causing the computer to execute the above-mentioned Gaussian process regression method for system identification.

[0060] The advantages and benefits of the present application are as follows:

[0061] ​1. The application obtains a data subset uniformly distributed without complex prior knowledge, reduces the size of the data subset under the premise of retaining data reflecting the key features of the model, thereby effectively reducing the calculation amount of Gaussian process regression;

[0062] 2. The application has strong adaptive ability. Since the data screening is based on the Euclidean distance, it is suitable for data of any dimension, and the algorithm can online calculate the volume of the space covered by the data, thereby adaptively adjusting the Gaussian kernel parameter γ to ensure good screening effect.

[0063] 3. The application can promote the use of Gaussian process regression in automatic control field. For the commonly used data with uneven spatial distribution in automatic control field, the application can extract a small-scale data subset with uniform distribution for the training of Gaussian process regression, thereby effectively reducing the calculation amount, expanding the use range of Gaussian process regression method in automatic control technology, and promoting the realization of more accurate and stable controller. BRIEF DESCRIPTION OF DRAWINGS

[0064] Figure 1 The figure is a whole flow chart of a Gaussian process regression method for system identification in an embodiment of the application.

[0065] Figure 2 The figure is a state motion curve of a control system in a specific embodiment of the application.

[0066] Figure 3 The figure is a distribution diagram of the input of the original data set in a two-dimensional plane in a specific embodiment of the application.

[0067] Figure 4 The figure is a distribution diagram of the extracted data subset in a specific embodiment of the application.

[0068] Figure 5 The figure is a distribution diagram of the test set in a specific embodiment of the application.

[0069] Figure 6 The figure is a distribution diagram of the prediction error and variance estimation value of the Gaussian process regression in a specific embodiment of the application

[0070] Figure 7 The figure is a frequency distribution diagram of the prediction error of the Gaussian process regression in a specific embodiment of the application DETAILED DESCRIPTION

[0071] The application proposes a Gaussian process regression method and device for system identification, which is described in further detail as follows in combination with the drawings and a specific embodiment:

[0072] The embodiment of the application proposes a Gaussian process regression method for system identification, which comprises:

[0073] The original data of the model to be identified in the automatic control system are obtained to form an original dataset. By assigning a corresponding importance factor to each set of original data, the importance factor is updated based on the Euclidean distance between the input data in the original data. The original data with low importance in the original dataset is removed to obtain a uniformly distributed subset of data.

[0074] The input data of the data subset is used to form a training input matrix, and the output data of the data subset is used to form a training output matrix. By performing Gaussian process regression training on the data subset, a covariance matrix considering data noise is obtained.

[0075] Based on the training output matrix, the training output matrix and the covariance matrix, the output prediction result corresponding to the input data to be predicted in the model to be identified is calculated to achieve system identification.

[0076] In a specific embodiment of the present invention, the overall process of the Gaussian process regression method for system identification is as follows: Figure 1 As shown, it includes the following steps:

[0077] 1) Obtain the original data of the model to be identified in the automatic control system to form the original dataset. Assign an importance factor to each group of original data in the original dataset, update the importance factor based on the Euclidean distance between the input data in the original dataset, and remove original data with low importance (i.e., those that are close in distance) to obtain a uniformly distributed subset of data. The specific steps are as follows:

[0078] 1-1) Obtain the original data of the model to be identified in the automatic control system to form the original dataset.

[0079] In this embodiment, the input dimension n is set according to the dimension of the model data to be identified. i Output dimension n o (This method is applicable to any input dimension n) i and output dimension n o For example, if the model to be identified is in the form of y = f(x1, x2), x1, x2, y ∈ R, then the input dimension n should be set. i =2, output dimension n o =1.

[0080] In this embodiment, the original dataset is denoted as S = {(x (j) ,y (j) ),j=1,2,...,P}, where and These are the input and output data for the j-th data set, respectively, and P is the total number of data sets. The input and output data in the original dataset are known, and there are no special requirements for the value of P in this embodiment.

[0081] 1-2) Create a data subset D initially empty, with a capacity C, which should not be too large or too small. If too large, the Gaussian process regression will have a large amount of computation; if too small, the extracted data subset will not be able to fully represent the system model, and the specific value can be obtained according to the model complexity and debugging.

[0082] 1-3) Let j = 1;

[0083] 1-4) Determine the number of stored data groups N in the data subset D:

[0084] If N = 0, it means that the current data subset D is empty, add the jth group of data (x (j) ,y (j) ) to D and set the importance factor μ j = 1, let N = N + 1, and record the maximum value M i of the ith element in the first j groups of input data in the original data set, where when j = 1, i = 1, 2, …, n i ; set the minimum value m i of the ith element in the first j groups of input data in the original data set, where when j = 1, i = 1, 2, …, n i ; for M i and m i , there is the following relationship: i = 1, 2, …, n i , k = 1, 2, …, j; take the currently added x (j) as the last added input data x last in the current data subset D, i.e. x last = x (j) , and then go to step 1-6).

[0085] If N > 0, it means that the current data subset D is not empty, update the maximum and minimum values of each element in the first j groups of input data in the original data set:

[0086]

[0087] Where “:” represents assignment.

[0088] Then calculate the Euclidean distance d between the input data x (j) in (x (j) , y (j) ) and the last added input data x last in the current data subset D:

[0089] d = ||x (j) -x last ||2 (2)

[0090] Let the set of the first j input data in the original data set be X = {x (k) |k = 1,2,...j}, calculate the minimum n i -rectangle (k-rectangle is defined as R k space that satisfies a i ≤x i ≤b i (1≤i≤k),a i ,b i ∈R, where a i is the minimum value of the ith input element and b i is the maximum value of the ith input element) volume V:

[0091]

[0092] Then calculate the receivable distance d a :

[0093]

[0094] Then go to step 1-5);

[0095] 1-5) Determine the Euclidean distance:

[0096] 1-5-1) If d≤d a , it means that the input data x (j) of this round is too close to the input data x last last added to the current data subset D, and is less important, so the jth original data is not added to the current data subset D.

[0097] Then determine j: if j = P, all data in the original data set has been traversed, and the current data subset D is the final data subset D, then go to step 2); if j

[0098] 1-5-2) If d>d a , add the jth original data (x (j) ,y (j) ) to the current data subset D and set the importance factor μ j of this data to 1, let N = N + 1, and take the current added x (j) as the input data x last last added to the data subset D, i.e. let x last=x (j) .

[0099] Then update the Gaussian kernel parameter γ using the following formula when V > 0:

[0100]

[0101] Where s is the standard deviation of the Gaussian kernel. When V = 0, we take γ = 0.

[0102] It should be noted that s can only be calculated if data exists in D.

[0103] Calculate the value of each data set (x) in the current subset D. (k) ,y (k) The importance of reducing gain G k k≠j:

[0104]

[0105] Update each data set (x) in the current subset D. (k) ,y (k) Importance factor μ k k≠j:

[0106] μ k :μ k ·G k (7)

[0107] Then proceed to steps 1-6);

[0108] 1-6) Determine j:

[0109] If j = P, then all data in the original dataset has been traversed, and the current data subset D is the final data subset D. Then proceed to step 2).

[0110] If j < P, then determine N: if N > C, it means the data subset is full, and the data with the smallest importance factor (x) in the current data subset D is selected. (r) ,y (r) ), Remove D, then let j = j + 1, and then return to step 1-4); if N ≤ C, then let j = j + 1, and then return to step 1-4).

[0111] 2) Perform Gaussian process regression training on the data subset D obtained in step 1), and use the training results to predict the output data corresponding to the input data to achieve system identification. The specific steps are as follows:

[0112] 2-1) Gaussian process regression training on the data subset D obtained in step 1); taking the final data subset D obtained in step 1) as the training set, the number of data groups of the training set (i.e. the number of data groups in the data subset D) is N, the training input matrix of the training set is X train = [x (1) ,x (2) ,...,x (N) ] T , and the training output matrix is Y and are column vectors, representing the i-th input data and output data in the data subset D respectively, i = 1, 2,..., N.

[0113] The covariance matrix K is calculated using the square exponential kernel y :

[0114]

[0115] where K y is the covariance matrix considering data noise, K train is the covariance matrix without considering data noise, is the noise variance of the output data, which should be set according to the actual noise size, I is an N x N identity matrix, and the square exponential kernel parameter p can be set according to the prior knowledge of model complexity and by debugging.

[0116] The matrices X train , Y train , K y are saved, which are needed for calculation during prediction.

[0117] 2-2) Prediction is performed using the results of step 2-1) to achieve system identification.

[0118] In this embodiment, it is assumed that n input data j = 1, 2,..., n corresponding output vectors need to be predicted. First, the covariance matrix K * between the training data and the prediction data is calculated, and the covariance matrix K ** between the prediction data is calculated:

[0119]

[0120] Then, the prediction matrix Y * and the covariance matrix Σ * of the output are calculated:

[0121]

[0122] where is the covariance matrix corresponding to the predicted value, while the covariance matrix * The jth element on the diagonal is the predicted value the variance estimate But the calculated value here is located in the interval [0, 1), which can only represent the relative size and has not been scaled to the real variance. Since the variance estimate is located in the interval [0, 1), it can be easily converted into a confidence level, thereby helping to improve the stability of the GPR prediction value when used for control. Assuming that the present application is applied to system identification and adaptive control, i.e., using Gaussian process regression to learn the system model and using the model prediction value and the predicted variance for the control law, through certain reasonable design, a high-precision and stable control effect can be achieved.

[0123] The method described in the present application will be further described in detail below in conjunction with a specific embodiment.

[0124] One specific embodiment of the present application adopts a digital simulation example in the field of automatic control:

[0125] In this example, for the control problem of a second-order system described by the following equation:

[0126]

[0127] where l, represent the displacement, velocity, and acceleration of the system, respectively, and u represents the control quantity. represents the dynamic model of the system, which has different forms for different actual systems. The form in equation (11) is only for the convenience of describing the embodiment of the present application.

[0128] In many control scenarios, the dynamic model cannot be accurately analytically modeled, and only a poor-performing model-free control method can be used. If the system identification can be performed using the motion data of the system, a better-performing model-based control method can be used.

[0129] In one specific embodiment of the present application, the Gaussian process regression method for system identification includes the following steps:

[0130] 1-1) Obtain the original data set:

[0131] In this embodiment, in order to identify the model the original data of the displacement l, the velocity the acceleration and the control quantity u need to be obtained. These data can be obtained through sensors or observers.

[0132] In order to obtain motion data with good spatial coverage, the displacement l of the system can track a target signal l in the form of random steps d In this embodiment, the step interval of the random step target signal is 5 s, and the signal size after each step is randomly selected from the interval [-5 m, 5 m]. Figure 2 The state motion curve of the control system in one specific embodiment of the present application. Figure 2 The state change of the system within 20 s is shown. Figure 2 The solid line in (a) represents the target signal l d over time.

[0133] Since is unknown, only a low-precision control using a model-free control method can be used to obtain motion data, and in this embodiment, a PID controller is used for control. The control law and the PID parameters obtained after debugging are as follows:

[0134]

[0135] wherein l d is the target signal of displacement, e is the tracking error of displacement, K p is the proportional gain, K i is the derivative gain, and K d is the integral gain.

[0136] wherein Figure 2 The solid line in (a) represents the target signal l d over time, and the dashed line represents the change of displacement l over time. The horizontal axis represents time, and the vertical axis represents displacement. Figure 2 (b) shows the speed of the system over time. The horizontal axis represents time, and the vertical axis represents speed. Figure 2 As can be seen from (a), the displacement of the system does not track the target signal well, because a model-free control method with poor performance is used.

[0137] Since the model to be identified is For the convenience of subsequent description, let the vector be the input, and the scalar be the output, so the model to be identified is In this embodiment, 200 s of input data and output data are collected to form an original data set, with a sampling frequency of 100 Hz, i.e. 20000 groups of data are collected, wherein white noise with a standard deviation of 0.01 is added to each output data.

[0138] The obtained original data set is represented as S={(x (j) , y(j) Let {j = 1, 2, ..., 20000}, where the superscript (j) represents the j-th data group. The original dataset input in this embodiment is distributed in a two-dimensional plane as follows: Figure 3 As shown in the figure, the horizontal axis represents the first component x1 of the input data, and the vertical axis represents the second component x2 of the input data. Each point in the figure represents one input data point. The figure demonstrates that the data has good spatial coverage and can be used for system identification.

[0139] Because the model to be identified in this embodiment It is a two-input, one-output function, so the input dimension n is set. i =2, Output dimension n o =1.

[0140] 1-2) Create an initially empty data subset D and set the size C of the data subset.

[0141] In this embodiment, C = 200; the total number of stored data N = 0.

[0142] 1-3) Let j = 1;

[0143] 1-4) Determine the number of data sets N already stored in the data subset D:

[0144] If N = 0, it means that the current data subset D is empty, and the j-th data (x (j) ,y (j) Add to D and set the importance factor μ for this set of data. j =1, let N = N+1, and denote the maximum value of the i-th element in the input data from the 1st to the j-th group in the original dataset as i = 1, 2, minimum value When i = 1, 2, the following relationship exists: i = 1, 2, k = 1, 2, ..., j; Add the currently added x (j) x is the last input data added to the current subset D. last , that is, x last =x (j) Then proceed to steps 1-6).

[0145] If N > 0, it means that the current data subset D is not empty. First, the maximum and minimum values ​​of each element in the input data from the 1st to the jth group of the original dataset are:

[0146]

[0147]

[0148] i = 1, 2

[0149] where ":" denotes assignment.

[0150] Then calculate the Euclidean distance d between the input x (j) in (x (j) ,y (j) ) and the last added input data x last in the current data subset D:

[0151] d = ||x (j) - x last ||2

[0152] Let the set consisting of the first j groups of input data in the original data set be X = {x (k) | k = 1, 2,... j}, calculate the minimum n i -dimensional hypercube in R k that contains X (k-hypercube is defined as R i that satisfies a i ≤ x i ≤ b i (1 ≤ i ≤ k), a i , b a ∈ R, such as 1-hypercube is a closed interval, 2-hypercube is a rectangle, and 3-hypercube is a cuboid), and the volume V of the hypercube:

[0153]

[0154] Then calculate the acceptable distance d a by the following empirical formula obtained by experiments:

[0155]

[0156] Then go to step 1-5);

[0157] 1-5) Determine the Euclidean distance:

[0158] 1-5-1) If d ≤ d (j) , it means that the input data x last of this round is too close to the last added input data x a in the current data subset D, and is less important, so the jth group of data is not added to the current data subset D;

[0159] Determine j: if j = P, all data in the original data set have been traversed, and the current data subset D is the final data subset D, then go to step 2); if j < P, let j = j + 1, and then return to step 1-4);

[0160] 1-5-2) If d > d a , add the jth group of data (x(j) , y (j) ) is added to D and the importance factor μ of this group of data is set j = 1, N = N + 1, and the currently added x (j) is recorded as the input data x last which was last added in the data subset D, i.e. let x last = x (j) .

[0161] Then the Gaussian kernel parameter γ is updated as follows:

[0162]

[0163] where s is the standard deviation of the Gaussian kernel. When V = 0, take γ = 0.

[0164] The importance reduction gain G (k) , y (k) ) of each group of data (x k , y (k) ) in the current data subset D is calculated, k ≠ j:

[0165]

[0166] The importance factor μ of each group of data (x (k) , y k ) in the current data subset D is updated, k ≠ j:

[0167] μ k : μ k · G k

[0168] Then go to step 1-6);

[0169] 1-6) Determine j:

[0170] If j = P, then all data in the original data set have been traversed, and the current data subset D is the final data subset D, then go to step 2);

[0171] If j < P, then determine N: if N > C, then the data subset is full, remove the data with the smallest importance factor from the current data subset D, and record this data as (x (r) , y (r) ), Then let j = j + 1, and then return to step 1-4); if N ≤ C, then let j = j + 1, and then return to step 1-4).

[0172] In this embodiment, in order to reflect the data subset extraction effect of the method described in this embodiment, the same size of data subset is extracted from the same original data set by random selection method as a control, and the effects of the two methods are shown in Figure 4 . Among them, Figure 4 (a) represents the data subset obtained by random selection method, Figure 4 (b) represents the data subset obtained by using the method of this embodiment Figure 4 (a) and Figure 4 (b) The abscissa represents x1 and the ordinate represents x2, and each point in the figure represents an input data. From Figure 4 (a), it can be seen that the data subset obtained by random selection method is mainly concentrated near x2=0, which retains the distribution characteristics of the original data set, and has poor space coverage, which is not conducive to system identification; from Figure 4 (b), it can be seen that the data subset obtained by using the method of this embodiment is evenly distributed, has good space coverage, and is conducive to system identification.

[0173] 2) Gaussian process regression training is performed on the data subset D obtained in step 1), and the output data corresponding to the input data is predicted by using the training result, so as to realize system identification

[0174] 2-1) Gaussian process regression training is performed on the data subset D obtained in step 1); the input matrix of the training data is denoted as X train =[x (1) ,x (2) ,...,x (N) ] T , and the output matrix is Y train =[y (1) ,y (2) ,...,y (N) ] T . The square exponential kernel is used to calculate the covariance matrix K y considering data noise:

[0175]

[0176]

[0177]

[0178] where K y is the covariance matrix considering data noise, K is the covariance matrix without considering data noise, is the noise variance of the output data, and in this embodiment, σ d=0.01, I is an N×N identity matrix, and the squared exponent kernel parameter ρ can be set according to prior knowledge of model complexity and through debugging. In this embodiment, it is set to 10.

[0179] Save matrix X train ,Y train ,K y It is needed for calculations during prediction.

[0180] 2-2) Use the results of step 2-1) to make predictions in order to achieve system identification.

[0181] In this embodiment, to verify the fitting effect of Gaussian process regression, a grid point set on the input plane can be taken as a test set, and the fitting effect can be tested by the prediction accuracy of Gaussian process regression on the test set. Figure 3 As can be seen from the data, the first component x1 of the training data input is located in the interval [-5, 5], and the second component x2 is located in the interval [-70, 90]. Therefore, in this embodiment, a grid set of 25 points at equal intervals in each dimension, distributed on a rectangle with -5≤x1≤5 and -70≤x2≤90, is selected as the test set.

[0182] Figure 5 This is a distribution diagram of the test set in a specific embodiment of the present invention, where the horizontal axis represents x1 and the vertical axis represents x2. Figure 5 Each point in the set represents a test point. Let the test set be denoted as . Let j represent the j-th test point, and the total number of test points is n = 625.

[0183] Next, the predicted values ​​of the Gaussian process regression on the test set are calculated. First, matrix K is calculated. * and K ** :

[0184]

[0185]

[0186] Then calculate the output prediction matrix Y. * Covariance Σ * :

[0187]

[0188]

[0189] in It corresponds to The predicted value, and the covariance matrix Σ * The j-th element on the diagonal is the predicted value. variance estimate However, the calculation obtained here... It lies on the interval [0,1), and can only represent relative size; it is not scaled to the true variance. Let y * and σ 2 When applied to control laws, it is beneficial for achieving high-precision and stable control. For example: for The predicted value y * It can be used for calculating control quantities, such as setting control quantities. The control system can then be approximated as a linear second-order system. The damping ratio ξ and natural frequency ω are set using classical control theory. n High-precision control can be achieved. For the variance estimation σ of the predicted value... 2 It can be used to judge the accuracy of the predicted value and to adjust the control law. For example, in areas not covered by the training data, the prediction error will be relatively large, and this can be addressed by adjusting the variance σ. 2 It is known that at this time y * No longer approximate A model-free control method should be adopted to ensure the stability of the system.

[0190] In contrast, this embodiment also uses a subset of data extracted by random selection for GPR training and prediction. Figure 6 This illustrates the prediction error |yy| of GPR in a specific embodiment of the present invention. * |and variance estimate σ 2 The distribution. Looking at the direction of the rows, Figure 6 (a) Figure 6 (b) Figure 6 (c) shows the results of GPR training using a randomly selected subset of data. Figure 6 (d) Figure 6 (e) Figure 6 (f) shows the results of GPR training using a subset of data extracted using the Data Manager. Looking at the columns... Figure 6 (a) Figure 6 (d) shows the distribution of the training data, and the image content is consistent with... Figure 4 The same is placed here for comparison; Figure 6 (b) Figure 6 (e) shows the distribution of GPR prediction errors. The horizontal axis of each plot represents x1, and the vertical axis represents x2. The color intensity of the contour lines indicates the magnitude of the prediction error; darker colors indicate smaller errors. This can be seen through comparison. Figure 6(e) The more dark color areas, the better the model fitting effect after GPR training using the data subset extracted by the method of the embodiment, and the light color areas with large prediction errors only appear in places where the training data does not cover; Figure 6 (c), Figure 6 (f) The figure shows the distribution of variance σ 2 , the horizontal coordinate of each figure represents x1, and the vertical coordinate represents x2. The color depth of the contour map represents the size of variance σ 2 . The darker the color, the smaller the variance. By comparing the images in the second column and the third column, it can be seen that there is a correlation between them, that is, the larger the variance, the larger the prediction error, which means that the size of the variance estimate value can represent the accuracy of the prediction. However, it is worth noting that variance is a statistical concept, which represents the possible size of the prediction error, so the image of the prediction error and the image of the variance are not one-to-one corresponding.

[0191] Figure 7 The figure shows the frequency distribution of the prediction error of GPR on the test set in one specific embodiment of the application, Figure 7 (a) is the result of GPR training using a randomly selected data subset, Figure 7 (b) is the result of GPR training using the data subset extracted by the method of the embodiment of the application. The horizontal coordinate of each figure represents the prediction error, and the vertical coordinate represents the frequency of the corresponding error segment, Figure 7 (a) The average value of the corresponding error is 25.333, Figure 7 (b) The average value of the corresponding error is 13.96. From Figure 7 It can be seen that the prediction error of GPR training using the data subset extracted by the method of the embodiment of the application is mainly concentrated below 10, while the prediction error of GPR training using a randomly selected data subset is generally large. Therefore, the method of the embodiment of the application can obtain better model recognition effect, thereby facilitating the application of the Gaussian process regression method in automatic control technology and realizing more accurate and stable control.

[0192] To achieve the above embodiment, the second embodiment of the application proposes a Gaussian process regression device for system identification, comprising:

[0193] A data subset construction module is configured to obtain a raw data set composed of raw data of a model to be identified in an automatic control system, update an importance factor corresponding to each raw data based on the Euclidean distance between input data in the raw data, eliminate raw data with low importance in the raw data set, and obtain a uniformly distributed data subset.

[0194] a Gaussian process regression training module, configured to form input data of the data subset into a training input matrix, form output data of the data subset into a training output matrix, and obtain a covariance matrix considering data noise by performing Gaussian process regression training on the data subset;

[0195] a prediction module, configured to calculate an output prediction result corresponding to to-be-predicted input data in the to-be-identified model based on the training output matrix, the training output matrix, and the covariance matrix, so as to realize system identification.

[0196] It should be noted that the foregoing embodiment of the Gaussian process regression method for system identification is also applicable to the Gaussian process regression device for system identification of the present embodiment, and will not be described herein again. According to the Gaussian process regression device for system identification provided in the embodiment of the present application, the original data set is obtained by acquiring original data of a to-be-identified model in an automatic control system, the original data with low importance is removed from the original data set by assigning a corresponding importance factor to each group of original data, updating the importance factor based on the Euclidean distance between input data in the original data, and obtaining a uniformly distributed data subset; the input data of the data subset is formed into a training input matrix, the output data of the data subset is formed into a training output matrix, and a covariance matrix considering data noise is obtained by performing Gaussian process regression training on the data subset; and the output prediction result corresponding to to-be-predicted input data in the to-be-identified model is calculated based on the training output matrix, the training output matrix, and the covariance matrix, so as to realize system identification. In this way, the uniformly distributed small sample set can be extracted from the commonly distributed uneven data in the automatic control field, thereby solving the problem of large calculation amount of Gaussian process regression when applied to system identification, and promoting the implementation of high-precision stable controllers.

[0197] To implement the foregoing embodiments, a third aspect of the present application provides an electronic device, comprising:

[0198] at least one processor; and a memory connected with the at least one processor in communication;

[0199] The memory stores instructions executable by the at least one processor, and the instructions are configured to perform the Gaussian process regression method for system identification.

[0200] To implement the foregoing embodiments, a fourth aspect of the present application provides a computer readable storage medium, which stores computer instructions for causing the computer to perform the Gaussian process regression method for system identification.

[0201] It should be noted that the computer-readable medium in the above disclosure can be a computer-readable signal medium or a computer-readable storage medium or any combination of the two. The computer-readable storage medium may, for example, but is not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, device or apparatus, or any combination of the above. More specific examples of a computer-readable storage medium can include, but are not limited to, an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disk read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the above. In the present disclosure, the computer-readable storage medium can be any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, device or apparatus. In the present disclosure, the computer-readable signal medium can include a data signal carried in a baseband or as a part of a carrier wave, which carries computer-readable program code. Such a propagated data signal can take many forms, including but not limited to an electromagnetic signal, an optical signal or any suitable combination of the above. The computer-readable signal medium can also be any computer-readable medium other than the computer-readable storage medium, which can send, propagate or transmit a program for use by or in conjunction with an instruction execution system, device or apparatus. The program code contained in the computer-readable medium can be transmitted by any suitable medium, including but not limited to a wire, a cable, an RF (radio frequency) or the like, or any suitable combination of the above.

[0202] The computer-readable medium described above can be contained in the electronic device described above; it can also exist separately and not be assembled into the electronic device. The computer-readable medium described above carries one or more programs, which, when executed by the electronic device, cause the electronic device to perform the Gaussian process regression method for system identification described above.

[0203] Computer program code for carrying out operations of the present disclosure can be written in any one or more programming languages or combinations of languages including object or visual programming languages such as Java, Smalltalk, C++ or conventional procedural programming languages such as the "C" programming language or similar programming languages. The program code can execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer can be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection can be made to an external computer (for example, through the Internet using an Internet Service Provider).

[0204] In the description of the specification, the description of the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" and the like means that the specific features, structures, materials or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present application. In the specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Also, the specific features, structures, materials or characteristics described can be combined in any appropriate manner in any one or more embodiments or examples. In addition, the person skilled in the art can combine and combine the different embodiments or examples described in the specification and the features of the different embodiments or examples, without contradiction.

[0205] In addition, the terms "first", "second", etc. are used only for the purpose of description and cannot be understood as indicating or implying relative importance or implicitly indicating the number of the technical features indicated. Therefore, the features defined with "first", "second" can explicitly or implicitly include at least one of the features. In the description of the present application, the meaning of "a plurality of" is at least two, for example, two, three, etc., unless otherwise explicitly specified.

[0206] Any process or method descriptions or descriptions of the flow diagrams in the specification or otherwise described herein can be understood as representing code modules, segments, or portions of code which include one or more executable instructions for implementing specific logic functions (or steps) in the process, and the various embodiments of the application can include additional or fewer functions. The various embodiments of the application can also modify other processes or methods described or otherwise suggested to perform the functions of the various embodiments of the application, and these modifications should be understood as falling within the scope of the various embodiments of the application. The various embodiments of the application should be understood to support a variety of additional modifications to otherwise result in performing the functions of the various embodiments of the application.

[0207] The logic and / or steps represented in the flowcharts and / or described herein, for example, can be considered as a sequence of instructions to implement logical functions, and can be embodied in any computer-readable medium for use by an instruction execution system, apparatus, or device, such as a computer-based system, processor- based system, or other system that can fetch the instructions from the instruction execution system, apparatus, or device and execute the instructions. For purposes of this specification, a "computer-readable medium" can be any apparatus that can contain, store, communicate, propagate, or transport the program for use by or in connection with the instruction execution system, apparatus, or device. The computer-readable medium can be a computer- readable storage medium or a computer-readable signal medium. The computer- readable storage medium can be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer-readable storage medium include the following: an electrical connection having one or more wires (electrical connections), a portable computer diskette (a magnetic device), a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, and a portable compact disc read-only memory (CDROM). In addition, the computer-readable medium can even be paper or another suitable medium upon which the program is printed, as the program can be electronically captured, for example, via optical scanning of the paper or other medium, then compiled, interpreted, or otherwise processed in a suitable manner, if necessary, and then stored in a computer memory.

[0208] It should be understood that aspects of the application can be implemented in hardware, software, firmware or combinations thereof. In the above embodiments, the various steps or methods can be implemented in software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented in hardware, as in another embodiment, any of the following technologies, known in the art, or their combinations can be used: discrete logic circuitry having logic gates for implementing logic functions on data signals, application specific integrated circuits having appropriate combinational logic gates, programmable gate arrays (PGA), field programmable gate arrays (FPGA), and the like.

[0209] Those skilled in the art can understand that all or part of the steps carried out by the above-mentioned embodiments can be completed by programs instructing related hardware, and the programs can be stored in a computer-readable storage medium. When the programs are executed, they include one or a combination of the steps of the method embodiments.

[0210] In addition, each of the functional units in the various embodiments of the present application can be integrated in one processing module, or each of the units can be physically present separately, or two or more units can be integrated in one module. The integrated module can be implemented in the form of hardware or in the form of a software functional module. When the integrated module is implemented in the form of a software functional module and sold or used as an independent product, it can also be stored in a computer readable storage medium.

[0211] The storage medium mentioned above can be a read-only memory, a magnetic disk or an optical disk, etc. Although the embodiments of the present application have been shown and described above, it should be understood that the above embodiments are exemplary and should not be construed as limiting the present application, and those skilled in the art can make changes, modifications, replacements and variations to the above embodiments within the scope of the present application.

Claims

1. A Gaussian process regression method for system identification, characterized in that, The method comprises the following steps: obtaining original data sets of a model to be identified in an automatic control system, updating importance factors of each set of original data based on Euclidean distances between input data in the original data, removing original data with low importance in the original data set, and obtaining a data subset with uniform distribution; composing a training input matrix by using input data of the data subset and composing a training output matrix by using output data of the data subset, and obtaining a covariance matrix considering data noise by performing Gaussian process regression training on the data subset; calculating output prediction results corresponding to input data to be predicted in the model to be identified based on the training input matrix, the training output matrix and the covariance matrix, and realizing system identification; wherein the obtaining of the data subset with uniform distribution comprises: 1) obtaining original data sets of a model to be identified in an automatic control system; Wherein, the input dimension of the original data of the model to be identified is , and the output dimension is ; Let the original data set be wherein and are the input data and the output data of the first group of data, respectively, is the total number of data groups. 2) Create a data subset initially empty with a capacity of ; 3) setting j=1; 4) on the data subset The number N of data sets stored in the memory is determined: If , it means that the current data subset is empty, the jth group of data is added to and the importance factor of the group of data is set to , let = +1, and let the maximum value of the ith element in the first to jth group of input data in the original data set be M i , where when j=1, ; let the minimum value of the ith element in the first to jth group of input data in the original data set be m i , where when j=1, ; ; the currently added is taken as the last added input data in the current data subset , i.e. , and then step 6) is entered; ; If , it means the current data subset is not empty set, update the maximum and minimum values of each element in the first group to the jth group of input data of the original data set as follows: （1） Computing the Euclidean distance between the input data in the current data subset and the last added input data in the current data subset : （2） Let the set of the first through jth input data sets of the original data set be , and calculate the volume of the smallest - hypercube in the space that contains . : （3） Computing receivable distance : （4） then entering step 5); 5) determining the Euclidean distance: 5-1) If then the jth group of raw data is not added to the current data subset ; then determining j: if j=P, all data in the original data set is traversed, and the current data subset D is the final data subset D; if j 5-2) If Then the j-th group of original data Add to current data subset The importance factor of this set of data is also defined. ,make = +1, adding the currently added As a subset of data The last input data added , that is to say ; The Gaussian kernel parameters are then updated as follows : When time: （5） wherein is the standard deviation of the Gaussian kernel; When time, take ; Compute current data subset Importance reduction gain for each group of data k≠j:​ （6） updating the current data subset the importance factor of each group of data k≠j:​ （7） then entering step 6); 6) determining j: if j=P, all data in the original data set is traversed, and the current data subset D is the final data subset D; If j < P, then determine N: if then it means that the data subset is full, and the data with the smallest importance factor in the current data subset D is removed and then j = j + 1, and then return to step 4); if N ≤ C, then j = j + 1, and then return to step 4).

2. The method of claim 1, wherein, the obtaining of the covariance matrix considering data noise comprises: The data subset As the training set, the number of data groups of the training set is , the training input matrix of the training set is , the training output matrix is , and are column vectors, respectively representing the i-th input data and the i-th output data in the data subset , i=1,2,…,N; Using a squared exponential kernel Computing the covariance matrix : （8） where is the covariance matrix considering data noise, is the covariance matrix not considering data noise, is the noise variance of the output data, is is the identity matrix, is the square exponential kernel parameter.

3. The method of claim 2, wherein, the calculating of the output prediction results corresponding to the input data to be predicted comprises: Let the input data to be predicted be , is the number of input data to be predicted; Computing a covariance matrix of training data and prediction data and prediction data : （9） The output prediction matrix is then calculated and covariance : （10） wherein , is a prediction value corresponding to , the covariance matrix the element on the diagonal of the th row and the th column is the variance estimate of the prediction value .

4. A Gaussian process regression apparatus for system identification, characterized by The method comprises the following steps: a data subset construction module is configured to obtain original data sets of a model to be identified in an automatic control system, update importance factors of each set of original data based on Euclidean distances between input data in the original data, remove original data with low importance in the original data set, and obtain a data subset with uniform distribution; a Gaussian process regression training module is configured to compose a training input matrix by using input data of the data subset and compose a training output matrix by using output data of the data subset, and obtain a covariance matrix considering data noise by performing Gaussian process regression training on the data subset; a prediction module is configured to calculate output prediction results corresponding to input data to be predicted in the model to be identified based on the training input matrix, the training output matrix and the covariance matrix, and realize system identification; wherein the obtaining of the data subset with uniform distribution comprises: 1) obtaining original data sets of a model to be identified in an automatic control system; Wherein, the input dimension of the original data of the model to be identified is , and the output dimension is ; Let the original data set be wherein and are the input data and the output data of the first group of data, respectively, is the total number of data groups. 2) Create a data subset initially empty with a capacity of ; 3) setting j=1; 4) on the data subset The number N of data sets stored in the memory is determined: like This indicates the current subset of data. It is an empty set, and the j-th group of data is... Add to The importance factor of this set of data is also defined. ,make = +1, let M be the maximum value of the i-th element in the input data from group 1 to group j in the original dataset. i When j=1, Set the minimum value of the i-th element in the input data from the first to the j-th group of the original dataset to be m. i When j=1, ; ; Add the currently added As a current subset of data The last input data added ,Right now Then proceed to step 6). If , it means the current data subset is not empty, update the maximum and minimum values of each element in the first group to the jth group of input data of the original data set as follows: （1） Computing the Euclidean distance between the input data in the current data subset and the last added input data in the current data subset : （2） Let the set of the first through jth groups of input data of the original data set be , compute the volume of the smallest - hypercube in that contains : （3） Computing receivable distance : （4） then entering step 5); 5) determining the Euclidean distance: 5-1) If then the jth group of raw data is not added to the current data subset ; then determining j: if j=P, all data in the original data set is traversed, and the current data subset D is the final data subset D; if j 5-2) If Then the j-th group of original data Add to current data subset The importance factor of this set of data is also defined. ,make = +1, adding the currently added As a subset of data The last input data added , that is to say ; The Gaussian kernel parameters are then updated as follows : When Time: （5） wherein is the standard deviation of the Gaussian kernel; When time, take ; Compute current data subset Importance reduction gain for each group of data Importance reduction gain for each group of data k≠j: （6） updating the current data subset the importance factor of each group of data the importance factor of each group of data k≠j: （7） then entering step 6); 6) determining j: If j = P, all data in the original data set is traversed, and the current data subset D is the final data subset D; If j < P, then determine N: if then it means that the data subset is full, and the data with the smallest importance factor in the current data subset D is removed and then j = j + 1, and then return to step 4); if N ≤ C, then j = j + 1, and then return to step 4).

5. An electronic device, comprising: Comprise: At least one processor; And the memory connected with the at least one processor in communication; Wherein, the memory stores instructions executable by the at least one processor, the instructions are set to execute the method of any one of claims 1-3.

6. A computer-readable storage medium, characterized in that, The computer readable storage medium stores computer instructions for causing the computer to execute the method of any one of claims 1-3.

Images