CN118296373A - A method for predicting cascading failures of infrastructure networks based on MPC nonlinear data-driven - Google Patents

Abstract

This paper discloses a method for predicting cascading failures of infrastructure networks based on MPC nonlinear data drive, which achieves full control of the system while balancing the input control data set and model control cost, focusing on spatial noise and discrete characteristics in real large-scale infrastructure networks. The effectiveness of the proposed scheme is verified using artificially synthesized networks and two real-world networks.

Description

Infrastructure network cascade failure prediction method based on MPC nonlinear data driving

Technical Field

The invention belongs to the technical field of the Internet of things, and particularly relates to an infrastructure network cascade failure prediction method based on MPC nonlinear data driving.

Background

The cascade failure is a key research hotspot in the field of complex networks, plays an important role in disaster early warning, disaster processing, fault recovery and the like, and is also one of main reasons for causing network paralysis of large-scale infrastructures. Simulation and control optimization of the network cascading failure process plays a significant role in controlling the propagation of cascading failures and in proposing a mitigation strategy for the cascading failures. However, the current research on the control of the cascade failure prediction model is mainly focused on the linear modeling of a network system, but the simple linear prediction model obviously cannot accurately reflect the relevant characteristics of a real network in consideration of the dynamic network characteristics and the node space propagation properties of a large infrastructure network. Therefore, a nonlinear model predictive control (Model Predictive Control, MPC) model oriented to cascade failure prediction is generated, and is more in line with the cascade failure propagation process in a real network environment.

Through summarizing the related researches in the past, the existing nonlinear data driving model is found to fully consider the dynamic characteristics of the network and the actual cascade failure propagation state, but noise and the discrete characteristics of a real system are often ignored when the dynamic characteristic information of the network is captured. And the existing controllable framework realizes complete control in theory, but does not solve the problem of reliable control and the problem of balance between the number of driving nodes and the control cost in a real large infrastructure network.

Therefore, how to realize cascade failure prediction in a large infrastructure network and balance model control reliability and cost is a difficult problem facing the safety of the large infrastructure network at present.

Disclosure of Invention

In order to overcome the defects of the prior art, the invention provides an infrastructure network cascade failure prediction method based on MPC nonlinear data driving, which realizes complete control of a system and balances input control data set and model control cost, and focuses on spatial noise and discrete characteristics in a real large infrastructure network. And uses a synthetic network and two real world networks to verify the validity of the proposed solution.

The technical scheme adopted for solving the technical problems is as follows:

Step 1: a multiple data collection stage;

step 1-1: controllable discrete linear system:

Wherein the system matrix The system state is thatThe system input isW (t) is system process noise, v (t) is system measurement noise, and y (t) is system output state;

assuming observable system state, the system matrix C is derived from the system interval matrix I and satisfies:

the interval matrix I specifies that each matrix element is bounded on the left and right sides by a left boundary I All possible values in between;

using zonotope model construction of the input matrix containing the presence of noise, assuming the zonotope constraint is satisfied for all system noise:

In the method, in the process of the invention, Zonotope, representing the location of noise w (t), which is defined by the center vectorAnd generating a matrixDefinition; Is a center vector describing noise, representing the center position of the noise; Is a generator matrix describing noise for defining the shape and size of the noise;

v (t) is expressed as:

In the method, in the process of the invention, Zonotope, representing where the noise v (t) is located, is defined by the center vectorAnd generating a matrixDefinition; Is a center vector describing noise, representing the center position of the noise;

Step 1-2: assume that the time step Av (t) satisfies the following constraint:

In the method, in the process of the invention, Zonotope, representing the time step Av (t), is represented by a center vectorAnd generating a matrixAnd (5) defining.Is a center vector describing the time step and represents the center position of the time step.Is a generator matrix describing a time step for defining the shape and size of the time step;

Consider K input-output trajectories of different lengths T _i, i=1 AndRepresenting sets of input and output states of the system under different trajectories, respectively, expressed as:

Wherein Y represents a system output state set, and U-represents a system input state set;

step 1-3: unfolding definition of Y:

Wherein Y ₊ and Y _- represent sets of boundary values for the output states of the system;

Total number of data points from all trajectories All available sets are represented by d= { U _-, Y }; for the system of formula (1), its input/output track d= { U _-, Y }, input state space matrixFurther described is:

In the method, in the process of the invention, A process noise matrix representing the system is used to describe the uncertainty of the system or the effect of noise on the output.A measurement noise matrix of the system for describing the effect of sensor errors or measurement noise on the output; a time step matrix representing the system for describing the effect of the time step on the output;

Step 2: a data optimization stage;

Step 2-1: modifying formula (1) to be:

y＝d_true(z_k)+w (9)

Using The regression function determines an unknown function s _true: Wherein the method comprises the steps of Is a gaussian noise which is a function of the noise,Is thatThe related characteristics needed to be used in regression testing;

Step 2-2: the prior s _true is modeled as a gaussian process:

Where m (z) and k (z, z') are any effective functions that calculate d _true means and covariance for a given Priori and training data dictionaryWith m observation pairs (z _i,y_i), there are:

The posterior distribution of the unknown function at test point z also follows gaussian distribution, with:

wherein:

using the zero mean function m (z) =0 and the average index kernel SEK gives:

Wherein the method comprises the steps of AndCovariance matrices of square sum and time step of output variances respectively;

the parameterized length is chosen to be positive half of the diagonal, resulting in:

meets the requirement that l _i is more than or equal to 0,

Thus, the free parameter of the core, i.e. the excess, is calculated as:

In the method, in the process of the invention, Representing a noise variance parameter in the kernel function for describing a noise level in the observed data;

step 2-3: to each dimension Combining together to obtain a complex made up of the given elementsThe approximation is:

Wherein the average value is And the variance Σ ^d (z) is calculated as:

Wherein, Representing the mean of the gaussian process at point z for each dimension separately,Representing the variance of the gaussian process at point z for each dimension separately;

step 3: model predictive control stage;

Step 3-1: let N be the length of the prediction boundary, AndAnd respectively obtaining input and output sequences of a system D= { U _- and Y }, and then performing convex quadratic optimization calculation to obtain a cost function:

Where Q _i = 0 and R _i＝0,…N-1 is a real symmetric positive definite matrix, Q _N is a real symmetric positive definite matrix, representing the weight matrix of the system final state y _N. The bias vectors for system state y _i and input u ₀, respectively. Further describing the objective as minimizing the distance between { y } _i and { y ^random}_i, the convex quadratic optimization calculation cost function reduces to:

where Q is a weight matrix that determines the importance of the observations in y; Representing a random state, which is a random observation compared with the system state y _i;

Order the The formula (20) is abbreviated as:

step 3-2: the MPC controller solves the following optimization problem at each time step of closed loop operation:

wherein ζ _c is the measured delay embedding vector and can be calculated as:

Wherein the method comprises the steps of And E _N define multidimensional states and input constraints, solve for optimal input sequences using Burgers' equationsThe problem translates into a locally optimal solution that solves an optimization problem of the form:

Wherein, F (x) is R ⁿ -R is an optimization objective function, and g (x) is R ⁿ→R^m is a constraint function;

Solving for an optimal solution for an input sequence To be obtainedThe method is applied to a complex network nonlinear system.

The beneficial effects of the invention are as follows:

The invention overcomes the defect that the traditional cascade failure model cannot balance the model control cost and the size of the input node data set, improves the cascade failure prediction accuracy and reduces the model control cost through a pure data-driven nonlinear prediction model.

Drawings

FIG. 1 is a flow chart of the overall process of the present invention.

FIG. 2 is a graph of the result of (a) minimizing the number of input edges in the ER network of the present invention; (b) results in the absence of noise.

FIG. 3 is a graph of the accuracy of the cascade failure prediction of the method of the present invention in four different propagation models, (a) experimental results in ER network; (b) results of experiments in US Powergrid networks; (c) results of experiments in OpenFlight networks.

FIG. 4 is a graph of the probability of state transition in four different propagation models for the method of the present invention, (a) experimental results in ER network; (b) results of experiments in US Powergrid networks; (c) results of experiments in OpenFlight networks.

Fig. 5 is a graph of the time cost of the present invention for different modes of peer-to-peer propagation in an ER network.

Detailed Description

The invention will be further described with reference to the drawings and examples.

The invention aims at realizing the balance problem of minimizing the input control data set and the model control cost on the premise of realizing complete control of a system aiming at the spatial noise and the discrete characteristic in a real large infrastructure network.

In order to achieve the above purpose, the technical scheme adopted by the invention is as follows:

Step 1: a multiple data collection stage;

Step 2: a data optimization stage;

step 3: model predictive control stage;

step 4: four control model instances.

Fig. 1 is an overall flow chart of the proposed framework of the invention.

(1) The specific implementation process of the step 1 is as follows:

because of the dynamic nature of real network environments, noise in the system has a great impact on model control, but previous research has generally used linearization of the system for simplicity of calculation, usually without taking system noise into account. Consider first a controllable discrete linear system:

Wherein the system matrix The system state is thatThe system input isAssuming observable system state, the system output matrix C is derived from the system interval matrix I and satisfies:

the interval matrix I specifies that each matrix element is bounded on the left and right sides by a left boundary I All possible values in between. Due to the presence of system noise, a multi-layer model is typically included that is consistent with the data [ A B ]. An input matrix containing the presence of noise is constructed using the zonotope model. Assuming that the limit of zonotope is satisfied for all system noise, then:

v (t) may be further expressed as:

further assume that all steps Av (t) satisfy the following constraint at all time steps:

the next goal is to solve the optimal control problem when the system model in equation (1) is unknown, but the input and noise output trajectories are available.

Consider K input-output trajectories of different lengths T _i, i=1And is also provided withI=1, …, K. All input data sequences can be expressed as:

further develop the definition of Y:

Total number of data points from all trajectories All available sets are denoted by d= { U _-, Y }. For a given system (equation (1)) input-output trajectory d= { U _-, Y }, the input state space matrix can be further described as:

this matrix contains all the matrices [ AB ] containing data and noise.

(2) The specific implementation process of the step 2 is as follows:

Modifying formula (1) to be:

y＝d_true(zk)+w (9)

determining an unknown function d _true using GP regression functions: Wherein the method comprises the steps of Is a gaussian noise which is a function of the noise,Is a relevant feature needed to be used in GP regression testing.

The a priori d _true is modeled as a gaussian process:

Where m (z) and k (z, z') are any effective functions that calculate d _true means and covariance for a given GP prior and training data dictionary And m observation pairs (z _i,y_i), there are:

The posterior distribution of the unknown function at test point z also follows gaussian distribution, with:

wherein:

here we give using a zero mean function m (z) =0 and an average index kernel (Square Exponential Kernel, SE):

Wherein the method comprises the steps of AndThe sum of squares of the output variances, respectively, the covariance matrix of the time step. We choose the parameterized length to be positive half-definite for the diagonal, resulting in:

meets the requirement that l _i is more than or equal to 0, Thus, the free parameter of the core, the so-called superstrate, can be calculated as:

So far, the GP regression process of one-dimensional output has been described. In general, one strategy to model GPs with multiple outputs is to treat each output dimension as an independent GP. Note that the super parameters are not necessarily equal and can be optimized independently. Thus, combining GP for each dimension results in a multiple GP approximation given as:

Wherein the average value is And variance Σ ^d (z) can be calculated as:

(3) The specific implementation process of the step 3 is as follows:

Next, in the MPC context, the control objective may be specified to minimize the cost function over a limited time range. A common approach is to use the model in equation (1) to predict the evolution of the system over a boundary and use these predictions to calculate the best input sequence for that boundary to minimize a given cost function. Then, only the first element of the calculated input sequence is applied to the actual system, thereby generating a new output value, and the whole process is repeated, which is sometimes referred to as the horizon fallback control technique. The method is described in detail below.

A significant feature of using the lifting linear predictor (equation (1)) is that although the original dynamics are nonlinear, the resulting MPC problem is a convex quadratic optimization problem (Quadratic Program, QP). Furthermore, the complexity of solving this problem if a so-called dense form is used can prove independent of lift magnitude, allowing the use of linear MPCs to efficiently solve this convex quadratic optimization problem.

Let N be the length of the prediction boundary,AndThe input and output sequence of the system input D= { U _-, Y } is the common convex quadratic optimization calculation cost function:

Where Q _i = 0 and R _i＝0,…N-1 is a real symmetric positive definite matrix. The cost function may be used to formulate a number of common control targets including tracking reference signals. For example, if the system is to be controlled, the target can be further described by { y ^random}_i by having its output measurements follow the indicated arbitrary time-dependent output sequence to minimize the distance between { y } _i and { y ^random}_i, so the cost function can be reduced to:

Where Q is a weight matrix that determines the importance of the observations in y. Note that in the above formula This term is not input-output dependent so it does not affect the optimal solution, i.e. it can be deleted directly in the representation. Order theEquation (20) can be abbreviated as:

This is also a special form of equation (19) and in the numerical example given by the present invention, mainly with this type of cost function, the MPC controller solves the following optimization problem at each time step of closed loop operation:

wherein ζ _c is the measured delay embedding vector and can be calculated as:

Wherein the method comprises the steps of And E _N define a multidimensional state and input constraints, which are typical convex quadratic optimization problems, the optimal input sequence can be solved using Burgers' equationThe problem translates into a (local) optimal solution that solves an optimization problem of the form.

Where f (x): R ⁿ → R is the optimization objective function and g (x): R ⁿ→R^m is the constraint function. Solving for an optimal solution for an input sequenceTo be obtainedThe method is applied to a complex network nonlinear system.

(4) The specific implementation process of the step 4 is as follows:

the proposed model is applied to four common signal dynamic propagation modes:

1) Random propagation: since the dynamic cascade failure process of the real network is irregular, the infectious disease model SIS, which is often used to simulate the nonlinear transmission process, is the most widely used. In this model, the infection corresponds to the failed node of the cascade failure, the susceptible population corresponds to the node that is affected first in the dynamic propagation of the cascade failure, Both forms are available. Because the dynamic process of the model is completely random, the probability of cascade failure is made to beAnd defining a cascade failure function alpha (l) as the probability that a failure-prone node is cascade-connected under the number of failed node neighbors l:

Pr(S→I∣l)＝α(l)＝1-(1-γ)^l (25)

Wherein gamma epsilon [0,1] is cascade failure propagation probability, that is, a node can be independently cascade connected to any failed node neighbor with gamma probability, and the recovery probability of cascade failure is defined as follows:

Pr(I→S)＝β (26)

the function of the cascade failure output node at this time is:

Wherein the method comprises the steps of Is the number of cascaded neighbor nodes to which node v _i is subject. Delta (x, y) is the Kronecker symbol (Kronecker delta). Note that for each case in equation (27), the result is a two-dimensional probability vector, where the first time is the probability that node v _i becomes/remains susceptible to cascading failure in the next time step, and the second term is the probability that the node becomes/remains in the failed state.

2) Complex propagation: the assumption of independent propagation in the SIS dynamics model is removed, and a complex cascade failure dynamics process is considered. The node result function is similar in form to equation (27), but the cascade failure function α (l) therein has a non-monotonic form, which can be re-expressed as:

where z (η) is a normalized failure node function such that α (l ^*)＝1.l^* and η >0 are parameters controlling the position of l.

3) Interaction propagation: cascade failure with interactions is defined as two interacting SIS kinetic models. In this case, there is s= [ S ₁S₂,I₁S₂,S₁I₂,I₁I₂ ] = {0,1,2,3}, similar to the simple SIS kinetic model, there isDefining a cascade failure probability function:

Wherein ζ.gtoreq.0 is a coupling constant, l _g is the number of neighbor nodes cascaded to by the failed node g, and a recovery probability of each cascade failure is defined as β _g, (g=1, 2). Where ζ >1 corresponds to the case where the failures are synergistic (i.e., after failure by one of the cascaded lines, the probability of being affected by the other cascaded line is increased), while if ζ <1 then contention is introduced (after failure has been affected by one of the cascaded lines, the probability of being affected by the other cascaded line is reduced). Case ζ=1 relies on two independent SIS dynamics that evolve simultaneously on the network. The result function consists of 16 entries, expressed as follows:

Where l _i,g is defined as the number of neighbor nodes of v _i cascade failure node g.

4) Community propagation: the model considers nodes to be composed of a fixed number of people N _i, including three states-a failure prone node (S), a failed node (I) or a restored node (R). Thus, the number of nodes in each state is tracked at a time. Further, the network G is weighted, and the weight represents the number of links between nodes. In this context,Is the average number of links from node v _j to node v _i. Finally, since the number of links between nodes is assumed to be stable, phi _i＝N_i is taken as a node attribute, and the node link proportion under dynamic transmission of each stream cascade failure is processed. More precisely, x _j＝(s_j,i_j,r_j) defines the state of node v _j, where s _j、i_j and r _j are the percentages of the failure-prone node, the failure node, and the restored node, respectively. From these definitions, the dynamic node result function is defined as:

wherein:

Where k _j is the time of node v _j. The function α (i, N) corresponds to the node failure probability iN each time step, i.e., any node iN the network is cascaded by neighbor nodes from the iN cascade failure nodes, and this probability can be expressed as:

Wherein R ₀ corresponds to the number of newly added failure nodes, and τ _r is the average recovery time of cascade failure.

Examples:

Experimental data set:

The software and hardware environment and the corresponding version numbers used in the invention are shown in table 1. The experimental data set used included three synthetic networks and 2 real large infrastructure networks (US Powergrid and Open Flight networks) with different topologies and presented their Degree, kcore and Triangle attribute values, respectively, with detailed network characteristics as shown in Table 2.

Table 1: experimental Environment configuration

Table 2 experimental data set-up

Experimental results:

(1) Influence of noise control on network modeling

According to the system model setting of formula (9), the output and input reference trajectories d= { U _-, Y }, the input and output constraintsAnd is also provided withWeight matrix R and Q zonotopes noiseAndModel cost functionModel observation delay ζ _c is used for respectively verifying optimal solutions of input sequences in four propagation modesIs effective in the following. Firstly, testing is carried out in an ER network of 2000 nodes, the situation that signals are randomly transmitted is eliminated, and the change situation that the number of required input edges is changed along with the transmission time of a cascade failure effect under the condition that whether noise exists or not is given, as shown in figure 2. Input set of three models in the presence of noise in ER networkIs about 5 times the maximum of the noiseless input set and community-based propagation is always kept low. In fig. 2, the community-based propagation method shows a larger value in the initial time step, because the community propagation is based on a higher modularity, and as the cascade failure continues to advance, the number of affected nodes increases, and the number of nodes required to be input for realizing the whole network model control also increases. When the time step t is less than 100, the node number required to be input for realizing model control is greatly reduced in three propagation modes, and the descending trend of interactive propagation is fastest. The three propagation modes achieve a substantially consistent slow increase in the number of nodes that the model control needs to input during the time step 100 to 200, at which point the number of remaining maximum communicating components in the network has smoothed. After the time step is longer than 200, the influence of cascade failure on the network is small, the network reaches a stable state, and the number of nodes which need to be input for realizing model control under three propagation modes is kept unchanged, but the number of edges which need to be added is integrally kept in a floating interval of [ -5,5] due to the influence of other noise factors of the network. The number of nodes required to be input for realizing model control under three different propagation modes after noise removal is shown in fig. 2 (b), the maximum control set number is about 1/5 of that under the existence of noise, and the performance of the proposed method is approximately consistent under the three propagation modes.

(2) Influence of propagation modes

Considering four common propagation modes of random propagation, complex propagation, interactive propagation and community propagation, the cascade failure prediction accuracy of the proposed model at the next moment in four different propagation models is respectively given in ER, US Powergrid and OpenFlight networks, and the test is developed by taking the system test time step as a reference, as shown in FIG. 3. The model shows higher prediction accuracy under four propagation modes along with the promotion of the cascade failure process, and particularly in a propagation mode based on communities. In the ER network in fig. 3 (a), the highest prediction accuracy can be achieved for all four propagation modes with a time step of about 460, and the node state in the network is constant at this time, that is, the network cascade failure will not affect the realization of the network functions and the normal operation of the network nodes. Because of the influence of the interactive network, as cascade failure advances, different networks are affected to different degrees, which also affects the accuracy of the proposed model cascade failure prediction, the cascade failure prediction accuracy based on interactive propagation shows one fluctuation in the time step 200, but maintains the same growing trend as other three propagation modes after 200 s. Because real large-scale infrastructure networks show high modularization characteristics such as water conservancy or electric power networks due to production and living requirements, a plurality of nodes always bear the main responsibility of the network, and the normal operation of the whole network is maintained. Once these critical nodes are attacked, the nodes connected to them will fail rapidly due to the highly dependent nature of the network nodes, and rapidly extend into the remaining nodes in the network, causing paralysis of the entire network. The performance of the proposed model is optimal in US Powergrid networks with a modularity of 0.94, and the prediction of cascading failure is fully achieved at t=400 s, as shown in fig. 3 (b). Because in a highly modular grid, the neighbors of the critical node at this time have been fully cascaded and failed, the set of nodes that need to be input for the entire network model control has been determined. The model presented in the OpenFlight network in fig. 3 (c) is substantially identical to the US Powergrid network, but the accuracy of the predictions for the four propagation modes is less than the results in the US Powergrid network.

The propagation state transition probabilities of the proposed model in four different propagation models are then given in the ER, US Powergrid and OpenFlight networks respectively, and the test is developed with reference to the system test time step, as shown in fig. 4. The cumulative distribution function (Cumulative Distribution Function, CDF), also called the distribution function, is the integral of the probability density function, and can fully describe the probability distribution of random variables in the proposed model, i.e. the set of edges that need to be input to achieve model control. A CDF map of the number of edges that need to be added in four propagation modes in three different networks is presented to show that the proposed model can minimize the set of model control inputs. First, in the ER network in fig. 4 (a), the same community-based propagation approach requires the least number of edges to be added and reaches steady state first, followed by interactive, complex and random based propagation. As shown in fig. 4 (b), in a highly modular US Powergrid network, the distribution of degrees is more uniform and reaches steady state earlier. The distribution state in the OpenFlight network of fig. 4 (c) is between the ER network and the US Powergrid network.

Finally, the time costs of the proposed model in three different propagation models are given in the ER, US Powergrid and OpenFlight networks respectively, and the test is developed with reference to the system test time step, as shown in fig. 5. The time consumption of the proposed model in the ER network in different propagation modes remains substantially the same. And the time step is more than 100 and then is stable, so that the scheme is proved to have lower time consumption while ensuring higher prediction accuracy for different cascade failure propagation modes of the network. The method has important significance for accurately protecting the normal nodes possibly affected by the cascade failure after the network is attacked, and saves the cost of predicting the cascade failure of the network. And time is saved for network failure recovery, and the probability of network cascade failure recovery is improved.

Claims (1)

1. An infrastructure network cascade failure prediction method based on MPC nonlinear data driving is characterized by comprising the following steps:

Step 1: a multiple data collection stage;

step 1-1: controllable discrete linear system:

Wherein the system matrix The system state is thatThe system input isW (t) is system process noise, v (t) is system measurement noise, and y (t) is system output state;

assuming observable system state, the system matrix C is derived from the system interval matrix I and satisfies:

the interval matrix I specifies that each matrix element is bounded on the left and right sides by a left boundary I All possible values in between;

using zonotope model construction of the input matrix containing the presence of noise, assuming the zonotope constraint is satisfied for all system noise:

In the method, in the process of the invention, Zonotope, representing the location of noise w (t), which is defined by the center vectorAnd generating a matrixDefinition; Is a center vector describing noise, representing the center position of the noise; Is a generator matrix describing noise for defining the shape and size of the noise;

v (t) is expressed as:

In the method, in the process of the invention, Zonotope, representing where the noise v (t) is located, is defined by the center vectorAnd generating a matrixDefinition; Is a center vector describing noise, representing the center position of the noise;

Step 1v2: assume that the time step Av (t) satisfies the following constraint:

In the method, in the process of the invention, Zonotope, representing the time step Av (t), is represented by a center vectorAnd generating a matrixDefinition; is a center vector describing a time step, representing the center position of the time step; is a generator matrix describing a time step for defining the shape and size of the time step;

Consider K input-output trajectories of different lengths T _i, i=1 AndRepresenting sets of input and output states of the system under different trajectories, respectively, expressed as:

Wherein Y represents a system output state set, and U _- represents a system input state set;

step 1-3: unfolding definition of Y:

Wherein Y ₊ and Y _- represent sets of boundary values for the output states of the system;

Total number of data points from all trajectories All available sets are represented by d= { U _-, Y }; for the system of formula (1), its input/output track d= { U _-, Y }, input state space matrixFurther described is:

In the method, in the process of the invention, A process noise matrix representing the system for describing the uncertainty of the system or the effect of noise on the output; a measurement noise matrix of the system for describing the effect of sensor errors or measurement noise on the output; a time step matrix representing the system for describing the effect of the time step on the output;

Step 2: a data optimization stage;

Step 2-1: modifying formula (1) to be:

y＝d_true(z_k)+w (9)

Using Regression function determines an unknown functionWherein the method comprises the steps ofIs a gaussian noise which is a function of the noise,Is thatThe related characteristics needed to be used in regression testing;

Step 2-2: the a priori d _true is modeled as a gaussian process:

where m (z) and k (z, z') are any effective functions that calculate d _true means and covariance for a given Priori and training data dictionaryWith m observation pairs (z _i,y_i), there are:

The posterior distribution of the unknown function at test point z also follows gaussian distribution, with:

wherein:

using the zero mean function m (z) =0 and the average index kernel SEK gives:

Wherein the method comprises the steps of AndCovariance matrices of square sum and time step of output variances respectively;

the parameterized length is chosen to be positive half of the diagonal, resulting in:

Satisfy the following requirements

Thus, the free parameter of the core, i.e. the excess, is calculated as:

In the method, in the process of the invention, Representing a noise variance parameter in the kernel function for describing a noise level in the observed data;

step 2-3: to each dimension Combining together to obtain a complex made up of the given elementsThe approximation is:

Wherein the average value is And variance Σ ^d (z) is calculated as:

Wherein, Representing the mean of the gaussian process at point z for each dimension separately,Representing the variance of the gaussian process at point z for each dimension separately;

step 3: model predictive control stage;

Step 3-1: let N be the length of the prediction boundary, AndAnd respectively obtaining input and output sequences of a system D= { U _- and Y }, and then performing convex quadratic optimization calculation to obtain a cost function:

where Q _i = 0 and R _i＝0,...N-1 is a real symmetric positive definite matrix, Q _N is a real symmetric positive definite matrix, representing a weight matrix of the system final state y _N; Offset vectors representing system state y _i and input u ₀, respectively; further describing the objective as minimizing the distance between { y } _i and { y ^random}_i, the convex quadratic optimization calculation cost function reduces to:

where Q is a weight matrix that determines the importance of the observations in y; Representing a random state, which is a random observation compared with the system state y _i;

Order the The formula (20) is abbreviated as:

step 3-2: the MPC controller solves the following optimization problem at each time step of closed loop operation:

wherein ζ _c is the measured delay embedding vector and can be calculated as:

Wherein the method comprises the steps of And E _N define multidimensional states and input constraints, solve for optimal input sequences using Burgers' equationsThe problem translates into a locally optimal solution that solves an optimization problem of the form:

Wherein f (x): r ⁿ -R is the optimization objective function, g (x): r ⁿ→R^m is a constraint function;

Solving for an optimal solution for an input sequence To be obtainedThe method is applied to a complex network nonlinear system.