Reachable set estimation method and system for periodic time-varying mechanical systems based on set intersection

Abstract

This invention relates to a method and system for estimating the reachability set of periodic time-varying mechanical systems based on set intersection. The method includes: establishing a state-space model of the periodic time-varying system; identifying and representing its periodic matrix as a sum of periodic time-varying terms; constructing a Lyapunov function based on linear piecewise division to obtain a time-varying criterion for reachability set estimation; constructing a multi-affine time-varying matrix polynomial; determining the convex constraints for reachability set estimation based on the properties of multi-affine, and designing an optimization function to obtain a constrained optimization problem to obtain the reachability set boundary; obtaining the optimal decay parameter through simulated annealing iterative optimization; selecting multiple decay parameters around the optimal decay parameter to obtain different boundaries, and generating a new reachability set boundary by taking the intersection of these parameters, and quantitatively calculating the degree of conservatism reduction of the new reachability set boundary. This invention can be used to estimate the reachability set of periodic time-varying mechanical systems, reduce the conservatism of reachability set estimation, and quantitatively evaluate the degree of conservatism reduction.

Description

Technical Field

[0001] The invention relates to a reachable set estimation method and system for a periodic time-varying mechanical system based on set intersection, and belongs to the field of periodic time-varying system control. Background Art

[0002] The reachable set of a dynamic system refers to the set that can constrain all system state trajectories under external interference conditions. The reachable set estimation of a time-varying system usually focuses on two goals: one is to find a finite boundary with certain mathematical rules. When the actual system is given an initial state and dynamic response characteristics, the trajectories of all time-varying state vectors are contained within the boundary; the other is to frame the system state within the smallest possible range to reduce the conservatism of the estimate and obtain a more accurate boundary while satisfying the first goal. For example, in the state space of a two-dimensional system, it is a prerequisite for the reachable set estimation that all states of the system are enclosed in the enclosed area, and measuring the size of the enclosed area enclosed by the reachable set determines the conservatism of the reachable set estimation and the accuracy of the estimated boundary.

[0003] In recent years, researchers have made outstanding contributions to the improvement of methods for obtaining reachable set constraints based on positive definite Lyapunov functions. For example, by constructing a co-positive Lyapunov function, the reachable set of a type of positive system with perturbation input is calculated. However, the application background of this method is mainly for discrete positive systems, and its conclusions are not applicable to continuous periodic time-varying systems. For another example, some scholars designed a state feedback controller for Markov switching systems based on the Lyapunov stability analysis principle and completed the reachable set estimation. From the simulation results of this method and its corresponding literature, it can be seen that this method cannot be directly applied to continuous periodic time-varying systems, and the corresponding Lyapunov attenuation parameters need to be manually selected, resulting in the estimated reachable set boundary being still conservative.

[0004] Considering that periodic phenomena and periodic characteristics are widely present in many engineering fields such as power systems, mechanical systems, network communications and control systems, it has become an urgent need to make the reachable set of periodic time-varying mechanical systems more accurately characterize the dynamic trajectory generated by the system in engineering applications and provide theoretical support for state tracking and verification of safety boundary reliability. Summary of the invention

[0005] The present invention provides a reachable set estimation method and system for a periodic time-varying mechanical system based on set intersection, aiming to solve at least one of the technical problems existing in the prior art.

[0006] The technical solution of the present invention relates to a reachable set estimation method and system for a periodic time-varying mechanical system based on set intersection. The method according to the present invention comprises the following steps:

[0007] A reachable set estimation method for a periodic time-varying mechanical system based on set intersection, characterized in that the method comprises the following steps:

[0008] S100, establishing a state space model of a periodic time-varying system, and expressing its period matrix in the form of sum of periodic time-varying terms through an identification method;

[0009] S200, constructing a Lyapunov function based on a linear piecewise form to obtain a reachable set estimation time-varying criterion based on the Lyapunov function, and then constructing multiple affine time-varying matrix polynomials;

[0010] S300, determining the convex constraint conditions of the reachable set estimation according to the multiple affine properties of the time-varying matrix polynomial, and obtaining the constrained optimization problem by designing the optimization function, and then solving the reachable set boundary;

[0011] S400, iteratively optimizing the attenuation parameter in the problem by a simulated annealing optimization method to obtain a numerical solution of the optimal attenuation parameter that can reduce conservatism;

[0012] S500, selecting multiple attenuation parameters around the optimal attenuation parameter to obtain different boundaries, generating a new boundary of the reachable set by taking the intersection method, and quantitatively calculating the degree of conservative reduction of the new boundary of the reachable set.

[0013] Further, in step S100,

[0014] T p is the basic period n x The state equation of the dimensional time-varying system is expressed as follows:

[0015]

[0016] In the formula, Represents n x dimensional state variables, is the periodic matrix obtained by identification, A(t)=A(t+T p ),Ψ(t)=Ψ(t+T p ), υ(t) is a bounded input signal, that is, there exists Make Established.

[0017] Further, in step S100,

[0018] T p The periodic time-varying parameter matrices A(t) and Ψ(t) are expressed as follows:

[0019]

[0020] In the formula, A0, Ah ,C h n x ×n x Order constant matrix, Ψ0,Ψ h ,Γ h n x ×n u Order constant matrix.

[0021] Further, the step S200 includes:

[0022] S201. Construct piecewise Lyapunov function V m (t) = x T (t)P m (t)x(t)(m=0,1,...,M-1), where P m (t) is a periodic symmetric positive definite matrix and adopts a linear piecewise form. p +mδ,lT p +(m+1)δ)P m (t) is expressed as:

[0023] P m (t) = P m +σ(t)P m+1

[0024] In the formula, the equal division parameter M represents the number of equal divisions in a cycle, and the width of the sub-interval is σ(t)∈[0,1) is a linear interpolation function;

[0025] S202, based on the known Lyapunov function reachable set estimation time-varying criterion, the non-convex time-varying function Q m (t)≤0 negative semidefinite problem is transformed into multiple affine time-varying matrix polynomials Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)] negative semidefinite problem.

[0026] Further, in step S202:

[0027] The known time-varying criterion Q for reaching set estimation based on Lyapunov function m (t) is expressed as follows:

[0028]

[0029] In the formula, About V m (t) is the Dini derivative; λ is the decay parameter of the Lyapunov function;

[0030] Among them, Qm The negative semidefinite form of the quadratic form of (t) is expressed as follows:

[0031]

[0032] Then the time-varying matrix polynomial [σ(t); α1(t),…,α H (t); β1(t),…,β H (t)] can be written as:

[0033]

[0034] In the formula, Λ 1,1,m ,…,Λ 1,2H+1,m and Λ 2,1,m ,…,Λ 2,2H,m They are all constant matrices.

[0035] Further, the step S300 includes:

[0036] S301, according to the properties of the affine function, a time-varying matrix polynomial is obtained. The necessary and sufficient condition for satisfying the following inequality is that the matrices corresponding to the endpoints of the time-varying parameter value range of the time-varying polynomial are all negative semidefinite;

[0037] Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)]≤0

[0038] According to the matrix polynomial Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)] is continuous and has a time-varying parameter σ(t)∈[0,1), and the polynomial limit set at σ(t)→1 is:

[0039]

[0040] S302, according to the lemma of time-varying reachable sets, P m (t) satisfies P m (t) = P m (t+T P )>0, and according to P m (t) Using linear equal interpolation, we can get P m The necessary and sufficient conditions for (t)>0 are:

[0041] P0>0,P1>0,…,P M-1 >0

[0042] S303, through the positive definite matrix Pm (t) is scaled, and the following formula is used as the optimization function:

[0043]

[0044] in, This is true for m=0,1,...,M-1;

[0045] According to the constraints Determine the reachable set boundary.

[0046] Further, the step S400 includes:

[0047] S401, using the cosine sampling function to construct the following attenuation parameters λ0,λ1,…,λ M-1 :

[0048]

[0049] Among them, the auxiliary parameter p∈[0,0.5] is set to make all the attenuation parameters λ m >0;

[0050] S402, obtaining a numerical solution of an optimal attenuation parameter λ that can reduce the conservatism of the reachable set by using a simulated annealing optimization algorithm;

[0051] Among them, the solution process of the simulated annealing algorithm is:

[0052] Initialize the decay parameter λ = λ (0) After that, the update law of the decay parameter is λ (k+1) =λ (k) +ρ·rand(0,1), where the maximum update step size ρ∈[0,1), rand(0,1) represents the random number generated in the interval [0,1]; when λ=λ (k+1) , substitute into the optimization function to calculate the optimal And calculate Δf= value; Set the current annealing temperature to T, if Δf<0 or Choose to accept λ=λ (k+1) , and set the next iteration temperature to μT, where the decay rate μ∈[0,1); repeat the above steps until the given convergence criterion is reached, and finally the optimal λ value can be obtained by numerical solution.

[0053] Further, the step S500 includes:

[0054] S501, select G different λ values ​​near the obtained optimal attenuation parameter λ and substitute them into the optimization problem to obtain G positive definite matrices Constructing new boundaries It is expressed as follows:

[0055]

[0056] Among them, the boundary equations L1, L2, ..., L G Respectively by get, is the result of the above different boundary equations L1, L2, ..., L G Find the new boundary obtained by intersection;

[0057] S502, according to the selected positive definite matrix Obtained by The determined ellipse area S1 or ellipsoid volume V1 is evenly scattered in the space to obtain the area S2 or volume V2 inside the reachable set, and then the efficiency of quantitatively obtaining the exact area or volume of the reachable set is obtained by the intersection method, which is expressed as follows:

[0058]

[0059] Where η S It represents the exact efficiency of the reachable set area relative to the original area after using the set intersection method; η V It represents the exact efficiency of the reachable set volume relative to the original volume.

[0060] The technical solution of the present invention also relates to a computer-readable storage medium on which program instructions are stored, and the above-mentioned method is implemented when the program instructions are executed by a processor.

[0061] The technical solution of the present invention also relates to a reachable set estimation system for a periodic time-varying mechanical system based on set intersection, the system comprising a computer device, a periodic time-varying mechanical system, and a system information identification and processing module and a plurality of sensors for feeding back system information to the computer device, the computer device comprising the above-mentioned computer-readable storage medium.

[0062] The beneficial effects of the present invention are as follows:

[0063] The present invention estimates the state reachable set of a periodic time-varying mechanical system based on the Lyapunov method, obtains a suitable attenuation parameter by searching and optimizing, so as to reduce the conservatism of the reachable set estimation result of the periodic time-varying mechanical system. According to the effectiveness of segmenting the periodic system based on the Lyapunov method, the present invention uses the Lyapunov matrix constructed based on the linear interpolation idea as the basis for analyzing the stability of the periodic time-varying system and obtaining the corresponding criterion. According to the properties of multiple affine scalar functions, the properties of affine functions can be extended from the scalar form to the characteristics of the time-varying matrix, solving the problem of strong coupling of multiple time-varying parameters in multiple affine polynomial constraints.

[0064] The present invention visualizes the reachable set estimation result of the periodic time-varying mechanical system as a plane geometric boundary generated by a two-dimensional piecewise curve or a solid geometric boundary generated by a three-dimensional point cloud. The generated boundary graphic area or volume size is compared with the ellipse (sphere) boundary graphic area or volume size before the set intersection, and a quantitative index is given to measure the reachable set estimation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 is a system structure diagram according to the method of the present invention.

[0066] Figure 2 is a flowchart of the method according to the present invention.

[0067] Figure 3 Schematic diagram of the system trajectory and reachable set boundary according to the first embodiment of the method of the present invention.

[0068] Figure 4 Schematic diagram of the system trajectory and reachable set boundary according to the second embodiment of the method of the present invention.

[0069] Figure 5 According to the second embodiment of the method of the present invention, a new three-dimensional reachable set boundary estimation result graph is obtained according to the sets corresponding to the boundaries of each ellipsoid. DETAILED DESCRIPTION

[0070] The concept, specific structure and technical effects of the present invention will be clearly and completely described below in combination with the embodiments and drawings to fully understand the purpose, scheme and effect of the present invention.

[0071] It should be noted that, unless otherwise specified, when a feature is referred to as being "connected" to another feature, it may be directly fixed or connected to another feature, or it may be indirectly fixed or connected to another feature. The singular forms "a", "said" and "the" used herein are also intended to include the plural forms, unless the context clearly indicates otherwise. In addition, unless otherwise defined, all technical and scientific terms used herein have the same meaning as those commonly understood by those skilled in the art. The terms used in this specification are intended only to describe specific embodiments and are not intended to limit the present invention. The term "and / or" used herein includes any combination of one or more related listed items.

[0072] The use of any and all examples, or exemplary language ("such as," "eg," etc.) provided herein, is intended merely to better illuminate embodiments of the invention and does not pose a limitation on the scope of the invention unless otherwise claimed.

[0073] Reference Figures 1 to 5In some embodiments, the reachable set estimation method for a periodic time-varying mechanical system with set intersection according to the present invention can reduce the conservatism of the estimation and improve the accuracy of the reachable boundary. The method includes at least the following steps:

[0074] S100, establishing a state space model of a periodic time-varying system, and expressing its period matrix in the form of sum of periodic time-varying terms through an identification method;

[0075] S200, constructing a Lyapunov function based on a linear piecewise form to obtain a reachable set estimation time-varying criterion based on the Lyapunov function, and then constructing multiple affine time-varying matrix polynomials;

[0076] S300, determining the convex constraint conditions of the reachable set estimation according to the multiple affine properties of the time-varying matrix polynomial, and obtaining the constrained optimization problem by designing the optimization function, and then solving the reachable set boundary;

[0077] S400, iteratively optimizing the attenuation parameter in the problem by a simulated annealing optimization method to obtain a numerical solution of the optimal attenuation parameter that can reduce conservatism;

[0078] S500, selecting multiple attenuation parameters around the optimal attenuation parameter to obtain different boundaries, generating a new boundary of the reachable set by taking the intersection method, and quantitatively calculating the degree of conservative reduction of the new boundary of the reachable set.

[0079] The present invention provides a method for estimating the reachable set of a periodic time-varying mechanical system, reducing the conservatism of the reachable set estimate, and quantitatively evaluating the degree of conservatism reduction, especially a method for visualizing the boundary and conservatism reduction of two-dimensional and three-dimensional reachable set estimates. The present invention constructs a mathematical model of a periodic time-varying system, designs a corresponding optimization function and linear matrix inequality constraints to obtain an optimization problem, and obtains the corresponding elliptical (spherical) boundary equation by solving the sub-optimization problem. After selecting different groups of parameters to obtain different boundaries, a new boundary equation is obtained by finding the intersection of the above boundary equations, and the degree of reduction in the conservatism of the reachable set can also be quantitatively approximated.

[0080] Specific implementation of step S100

[0081] The present invention establishes a state space model of a periodic time-varying system based on an identifiable periodic time-varying function, and analyzes and estimates the reachable set of the system through the Lyapunov method. It should be noted that in the process of modeling a periodic time-varying system, it is more common to decompose the periodic matrix in the state equation into a series of triangular series in the form of Fourier series expansion, that is, the coefficients are all in the interval [-1,1].

[0082] The present invention adopts the Fourier series method to model the periodic time-varying matrix, and expresses the periodic matrix in the form of adding periodic time-varying terms through the identification method. Among them, the periodic time-varying terms can be modeled by a series of trigonometric functions based on the Fourier method. Assume T p The time-varying continuous state equation with the basic period is:

[0083]

[0084] In the formula, Represents n x dimensional state variables, x(t) represents a series of common physical quantities, which can be used to represent a set of variables of the internal state of the system in state space theory. For example, in mechanical systems, basic physical quantities such as displacement, velocity, angular displacement, acceleration, etc. are usually selected as state variables. is the periodic matrix obtained by identification, A(t)=A(t+T p ),Ψ(t)=Ψ(t+T p ). The input signal υ(t) can represent physical quantities such as displacement and velocity, and the input signal is bounded, that is, there is a positive number Make Here, we take the example of identifying A(t) and Ψ(t) as a superposition of trigonometric function sequences composed of Fourier series, and obtain p The periodic time-varying parameter matrices A(t) and Ψ(t) satisfy the following form:

[0085]

[0086] Among them, A0,A h ,C h Yes x ×n x Order constant matrix, Ψ0,Ψ h ,Γ h Yes x ×n u The order constant matrices can be obtained based on the periodic time-varying system parameter identification or approximation. H (t) and β1(t),β2(t),…,β H (t) is in the interval [-1,1] with T p A continuous periodic time-varying function with a basic period.

[0087] It should be noted that the above modeling method of the present invention has a certain universality. It can not only determine the size of the time-varying function term H according to the actual scenario, but also the form of the periodic time-varying function can be more flexible. For example, the present invention adopts a modeling method based on Fourier series, which can be set in,

[0088]

[0089] In the formula, A pq (t) is the time-varying function of the system matrix A(t) with p rows and q columns. According to Fourier theory, the constant term is:

[0090]

[0091] The other higher-order harmonic factors a p,q,h ,c p,q,h Then record it as:

[0092]

[0093] Following the above process, after identification or approximation by Fourier method, the constant matrix ψ in ψ(t) can be obtained: h With Γ h The value of .

[0094] It is understandable that if α h (t), β h (t) The parameters are modified into other function forms that meet the conditions. The modeling method of the present invention is not limited to the form of Fourier series modeling.

[0095] Specific implementation of step S200

[0096] First, after the state space model of the periodic time-varying mechanical system based on the real object is determined, the piecewise Lyapunov function is constructed According to the given value of the equally divided parameter M, a basic period T p Divide into M equal intervals, where P m (t) is a periodic and symmetric positive definite matrix, that is, it satisfies P m (t) = P m (t+T p )>0. Specifically, since there are multiple coupled time-varying function terms in the time-varying constraint inequality in the reachable set estimation, the piecewise Lyapunov function V can be further determined after the number of segments and the attenuation parameter are determined. m (t), and construct a time-varying matrix polynomial with multiple affine properties based on the identification model of periodic systems and the time-varying reachable set lemma.

[0097] First, we use the periodic positive definite matrix P in linear piecewise form m (t)(m=0,1,...,M-1), and the equal division parameter M represents the number of equal divisions in a cycle, then the sub-interval width is Define σ(t)∈[0,1) as a linear interpolation function.p +mδ,lT p +(m+1)δ)

[0098]

[0099] So, the linear interpolation function P m (t) can be written as follows:

[0100] P m (t) = P m +σ(t)P m+1

[0101] In determining the time-varying function P m (t), the reachable set of a periodic time-varying mechanical system can be defined as:

[0102]

[0103] Then, according to the known time-varying criterion for reachable set estimation based on Lyapunov function, Q m (t)≤0, this non-convex time-varying function negative semidefinite problem is transformed into a multi-affine time-varying matrix polynomial Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)] negative semidefinite problem, and thus the state equation constraints of the reachable set of the periodic time-varying system are obtained by utilizing the properties of affine functions, which can improve the flexibility of the solution and facilitate the implementation of computer programs.

[0104] Specifically, if for every m=0,1,…,M-1 there exists a positive definite matrix P m (t) and λ m > 0, so that the matrix polynomial Q m (t)Satisfy:

[0105]

[0106] In the formula, It's about V m The Dini derivative of (t) is defined as:

[0107]

[0108] And at non-switching time, Substituting the above assumptions about the periodic time-varying matrices A(t) and Ψ(t), and combining them with the Lyapunov function V m Assuming (t), Q m (t) can be written as the following quadratic negative semidefinite form:

[0109]

[0110] And according to the time-varying matrix polynomial Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)] is formed by the state equation of the periodic time-varying mechanical system and the Lyapunov equation V on the mth subinterval m (t) is determined in its specific form, the time-varying matrix polynomial Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)] can be written as:

[0111]

[0112] In the formula, Λ 0,m , Λ 1,1,m ,Λ 1,2,m ,…,Λ 1,2H+1,m and Λ 2,1,m ,Λ 2,2,m ,…,Λ 2,2H,m are all constant matrices, and their corresponding constant matrices are:

[0113]

[0114]

[0115]

[0116]

[0117]

[0118]

[0119] In the formula, And the “*” represents the constant matrix that is omitted due to symmetry.

[0120] Specific implementation of step S300

[0121] After constructing the inequality conditions for time-varying polynomial decoupling, the convex constraints of the reachable set estimation can be determined according to the multiple affine properties of the matrix polynomial, and the conservatism of the reachable set estimation can be reduced by designing an optimization function. Specifically, the time-varying non-convex constraints in the original conditions are transformed into constant vertex constraints through the affine properties of the time-varying matrix, and then the optimization equation of the reachable set estimation can be obtained, and the elliptic (spherical) equation can be solved by the optimizer. Among them, the multiple affine properties are used to reveal the time-varying matrix Λ m [σ(t);α1(t),…,αH (t); β1(t),…,β H (t)] is related to the matrix values ​​of each endpoint in its domain.

[0122] First, according to the properties of affine functions, we know that Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)] is a time-varying matrix polynomial with multiple affine properties. Therefore,

[0123] Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)]≤0

[0124] The necessary and sufficient condition is that the matrices corresponding to the endpoints of the time-varying parameter range of the time-varying polynomial are all negative semidefinite. Among them, according to the fact that the time-varying matrix polynomial is continuous and σ(t)∈[0,1), the polynomial limit defined when σ(t)→1 is:

[0125]

[0126] Considering P m (t) is a periodic positive definite function and is in the form of linear equal interpolation when m = 0, 1, ..., M-1, we can get P m The necessary and sufficient conditions for (t)>0 are:

[0127] P0>0,P1>0,…,P M-1 >0

[0128] The present invention uses a positive definite matrix P m (t) Make appropriate scaling, i.e. And by P m (t) The determined reachable set closed curve is expressed as a boundary with certain geometric meaning. Specifically, if we directly use the time-varying positive definite matrix P m (t), then x T (t)P m (t)x(t)≤1 The boundary constraints on x(t) will contain time-varying terms and are difficult to solve. If is a constant positive definite matrix, then represents the boundary of the ellipse (ellipsoid) about x(t), so by constructing and meet according to It can be obtained by P m (t) must be a reachable set consisting of Constitute a subset of the reachable set. Among them, due to is a convex function (a convex constraint), then the reachable set estimation problem can be expressed as follows, with linear matrix inequality constraints and can be solved by convex optimization:

[0129]

[0130] Combined with the above derivation, we can see that due to the time-varying function Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)] and P m (t) are all piecewise matrix affine functions, so the condition for the negative semidefiniteness (positive definiteness) of the above time-varying matrix polynomial is the negative semidefiniteness (positive definiteness) at its endpoints. After constructing the constraints of the optimization problem, the present invention adopts the following formula as the optimization objective function:

[0131]

[0132] in It holds true for m=0,1,...,M-1, and is a symmetric positive definite matrix. On the one hand, this assumption can ensure that The generated reachable set must contain m (t) The generated reachable set; on the other hand, this hypothesis maximizes The determinant of can reduce the estimated area (volume) of the reachable set as much as possible.

[0133] Specific implementation of step S400

[0134] The present invention realizes parameter search through an optimization algorithm to obtain the optimal attenuation parameter value with reduced conservatism, thereby obtaining a new smooth boundary in the form of "set intersection" of reachable sets. It further reduces the conservatism of the reachable set boundary estimation by combining the optimization algorithm with the reachable set boundary equation to find the intersection, and approximates the range of reduced conservatism in a quantitative form.

[0135] Specifically, the cosine sampling function is first used to construct the following attenuation parameters λ0,λ1,…,λ M-1 , to reduce optimization variables and simplify the optimization problem:

[0136]

[0137] Among them, λ is the attenuation parameter of the Lyapunov function, and the auxiliary parameter p∈[0,0.5] is set to ensure that all the attenuation parameters λ m >0(m=0,1,...,M-1) holds, which means that the attenuation parameter is a positive number.

[0138] Then, the numerical solution of the optimal decay parameter λ that can reduce the conservatism of the reachable set is obtained through the simulated annealing algorithm. The simulated annealing algorithm solution process is as follows: Initialize the decay parameter λ = λ (0) After that, the update law of the decay parameter is λ (k+1) =λ (k) +ρ·rand(0,1), assuming that the maximum update step is ρ∈[0,1), rand(0,1) represents the random number generated in the interval [0,1]. (k+1) , substitute into the optimization function to calculate the optimal And calculate The value of. Let the current annealing temperature be T] if Δf<0 or Choose to accept λ=λ (k+1) The next iteration temperature is set to μT, where the decay rate μ∈[0,1). Repeat the above steps until the given convergence criterion is reached, and finally the optimal λ value can be obtained by numerical solution.

[0139] Specific implementation of step S500

[0140] The present invention adopts a metric optimization algorithm and a reachable set "boundary set intersection" form to obtain the degree of reduction of the reachable set, and quantitatively calculates the degree of reduction of the reachable set estimation conservatism through the degree of reduction of the reachable set.

[0141] Specifically, after obtaining the optimal value of the optimal attenuation parameter λ, G different λ values ​​are selected near the optimal value, and the selected G positive λ values ​​are substituted into the optimization problem for solution. The obtained symmetric positive definite matrix Recorded as Depend on The different boundary equations are L1, L2, ..., L G , thus obtaining a new boundary

[0142]

[0143] in, is composed of different boundary equations L1, L2, ..., L G The new boundary obtained by "stitching" is composed of the new boundary The enclosed area must be included in L1, L2, ..., L G Within the determined area, the new boundary of the reachable set based on the intersection of the set boundaries can be formed by intersecting multiple ellipses (spheres) obtained by intersection operations, thereby reducing the conservatism of the reachable set estimation and improving the accuracy of the reachable set boundary estimation.

[0144] Then, according to the new boundary L1, L2, …, L GThe smooth boundary obtained by "stitching" and satisfying Lipschitz continuity is available, and the area (volume) of the closed region enclosed by the new boundary can also be quantitatively estimated by the method of uniformly scattering points in space, thereby realizing quantitative calculation of the degree of conservative reduction in evaluating the reachable set.

[0145] Specifically, the above optimization problem is solved iteratively by the solver to obtain the symmetric positive definite matrix If the selected positive definite matrix is Then the symmetric positive definite matrix The area (volume) of the determined ellipse (ellipsoid) is recorded as:

[0146]

[0147] Then, the following scatter point method can be used to estimate the boundary The area (volume) of the bounding region:

[0148] Initialize a space W containing all system states, and evenly scatter points on the space to obtain the area (volume) inside the reachable set as S2(V2). After the boundary set intersection, perform the following calculations to quantitatively estimate the efficiency of estimating the reachable set area (volume):

[0149]

[0150] Where η S It represents the exact efficiency of the reachable set area relative to the original area after using the set intersection method; η V It represents the exact efficiency of the reachable set volume relative to the original volume.

[0151] Therefore, the present invention quantitatively calculates the degree of reduction in the conservatism of the reachable set estimate on the basis of reducing the coverage area of ​​the reachable set, and the size of the conservative reduction can be approximately quantitatively calculated by the interval scattering method.

[0152] Practical application test

[0153] The present invention conducts actual testing on the reachable set estimation system of periodic time-varying mechanical system based on set intersection. Figure 1 The method of the present invention can be applied to a periodic time-varying mechanical system. The information of the periodic time-varying mechanical system is obtained through a system information identification processing module and multiple sensors and fed back to a computer device, so as to obtain the reachable set estimation boundary of the above system by using the reachable set estimation method of the periodic time-varying mechanical system based on set intersection of the present invention. The reachable set estimation system of the present invention can not only output the reachable set estimation result of the system, but also quantitatively analyze the conservatism of the reachable set estimation and quantitatively calculate the degree of reduction of the conservatism of the reachable set.

[0154] In some embodiments, taking a second-order system as an example, the method of the present invention is applied to a second-order single-meshing gear system, where the fundamental period of the system is T p =2, the bounded disturbance is After the system is parameterized in state space, a periodic time-varying dynamic parameter matrix of a second-order system in the form of the present invention can be obtained. The equally divided parameter M=6 is taken, and three sets of parameters are used to obtain different elliptical boundaries, and the values ​​of the three sets of parameters are λ=1, p=0.25; λ=1.5, p=0.25 and λ=2, p=0.25.

[0155] See the reachable set estimation results for the second-order system example for Figure 3 , which shows the process of generating elliptical boundaries and reducing conservatism under three sets of parameter settings. Specifically, from the simulated annealing algorithm, it can be seen that when the auxiliary parameter p = 0.25 and the equal division parameter M = 6, the optimal value of the attenuation parameter λ obtained by iterative solution is λ = 1.402, so λ = 1, λ = 1.5, and λ = 2 are selected to generate different elliptical reachable set boundary estimation results, and the estimated system reachable set boundary is redrawn after set intersection calculation. It can be seen that the figure consists of 32 groups that meet The periodic time-varying trajectories generated by different perturbation inputs are all located inside all ellipse boundaries. The yellow block area represents the conservative area reduced by the parameter search optimization and set intersection estimation to obtain the new reachable set boundary compared to the ellipse boundary determined by λ=1.

[0156] Therefore, the present invention obtains appropriate attenuation parameters by searching and optimizing, and can obtain a series of elliptical boundary constraints with a small area and completely surrounding the system state when solving the two-dimensional reachable set optimization problem. By taking the intersection of each elliptical boundary, the conservatism of the reachable set estimation result of the periodic time-varying mechanical system can be reduced.

[0157] In some embodiments, the present invention is also applicable to a third-order system. The method of the present invention is applied to an electro-hydraulic servo vibration table system. The system model can be modeled as a third-order state equation and stabilized by state feedback. Based on the working conditions of the vibration table, according to the situation that the system parameter matrix is ​​periodic and time-varying, and there are certain bounded disturbances between the system state variables, the basic period T p =2, bounded disturbance After the system is parameterized in state space, a periodic time-varying dynamic parameter matrix of a third-order system conforming to the form described in the present invention is obtained, and then the equally divided parameter M=6 is taken, and four sets of parameters are used to obtain different ellipsoid boundaries, and the values ​​of the four sets of parameters are λ=1, p=0.25; λ=1.5, p=0.25; λ=0.5, p=0 and λ=0.5, p=0.25 respectively.

[0158] See the dynamic trajectory of the periodic time-varying system and the boundary estimation results of the ellipsoid reachable set based on the third-order system example for details. Figure 4 and Figure 5 , which shows the process of generating ellipsoid boundaries under four sets of parameter settings. Specifically, the interior of the point cloud contains periodic time-varying trajectories generated by 32 sets of different time-varying perturbations (but satisfying the same perturbation constraints), while the blue, green, red, and yellow point cloud boundaries correspond to the ellipsoid boundaries obtained by the four sets of parameters in the simulation example. Figure 5 Yes Figure 4 The sets corresponding to the boundaries of each ellipsoid take the common part of their point clouds for calculation, and obtain the new three-dimensional reachable set boundary estimation results. It can be seen that the periodic time-varying states generated by the system under the action of bounded disturbance are all contained in the interior of the reachable set point cloud boundary.

[0159] Therefore, the three-dimensional state reachable set estimation result obtained by the method of the present invention also reduces the conservatism compared with the traditional ellipsoidal reachable set, and the reduction in conservatism can be quantitatively calculated by reducing the volume.

[0160] From the above two example results ( Figures 3 to 5 ) It can be seen that the reachable set estimation result obtained by parameter search and set intersection processing in the present invention is compared with the reachable set in the form of an ellipse (sphere) obtained by directly estimating the given attenuation parameters. The area or volume covered by the boundary graphic is reduced to a certain extent, thereby improving the estimation accuracy of the reachable set boundary and reducing the conservatism of the reachable set estimation.

[0161] It should be appreciated that the method steps in the embodiments of the present invention can be implemented or implemented by computer hardware, a combination of hardware and software, or by computer instructions stored in a non-transitory computer readable memory. The method can use standard programming techniques. Each program can be implemented in a high-level process or object-oriented programming language to communicate with a computer system. However, if necessary, the program can be implemented in an assembly or machine language. In any case, the language can be a compiled or interpreted language. In addition, the program can be run on a programmed ASIC for this purpose.

[0162] Furthermore, the operations of the processes described herein may be performed in any suitable order unless otherwise indicated herein or otherwise clearly contradicted by context. The processes described herein (or variations and / or combinations thereof) may be performed under the control of one or more computer systems configured with executable instructions, and may be implemented as code (e.g., executable instructions, one or more computer programs, or one or more applications) that is executed collectively on one or more processors, by hardware, or a combination thereof. The computer program includes a plurality of instructions that may be executed by one or more processors.

[0163] Further, the method can be implemented in any type of computing platform that is operably connected to a suitable computer, including but not limited to a personal computer, a minicomputer, a mainframe, a workstation, a network or distributed computing environment, a separate or integrated computer platform, or in communication with a charged particle tool or other imaging device, etc. Various aspects of the present invention can be implemented in machine-readable code stored on a non-transitory storage medium or device, whether removable or integrated into a computing platform, such as a hard disk, an optical read and / or write storage medium, an RSM, a ROM, etc., so that it can be read by a programmable computer, and when the storage medium or device is read by the computer, it can be used to configure and operate the computer to perform the process described herein. In addition, the machine-readable code, or portions thereof, can be transmitted via a wired or wireless network. When such media includes instructions or programs that implement the steps described above in conjunction with a microprocessor or other data processor, the invention described herein includes these and other different types of non-transitory computer-readable storage media. When programmed according to the methods and techniques of the present invention, the present invention can also include the computer itself.

[0164] The computer program can be applied to input data to perform the functions described herein, thereby converting the input data to generate output data stored in a non-volatile memory. The output information can also be applied to one or more output devices such as a display. In a preferred embodiment of the present invention, the converted data represents physical and tangible objects, including specific visual depictions of physical and tangible objects produced on the display.

[0165] The above is only a preferred embodiment of the present invention. The present invention is not limited to the above implementation. As long as the technical effect of the present invention is achieved by the same means, any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the scope of protection of the present invention. Within the scope of protection of the present invention, its technical scheme and / or implementation method may have various modifications and changes.

Claims

1. A reachable set estimation method for periodic time-varying mechanical systems based on set intersection, characterized in that: The method is used to reduce the conservatism of the estimation and improve the accuracy of the reachable boundary, and the method comprises the following steps: S100, establishing a state space model of a periodic time-varying system, and expressing its period matrix in the form of sum of periodic time-varying terms through an identification method; S200, constructing a Lyapunov function based on a linear piecewise form to obtain a reachable set estimation time-varying criterion based on the Lyapunov function, and then constructing multiple affine time-varying matrix polynomials; S300, determining the convex constraint conditions of the reachable set estimation according to the multiple affine properties of the time-varying matrix polynomial, obtaining the constrained optimization problem by designing the optimization function, and then solving it to obtain the reachable set boundary; S400, iteratively optimizing the attenuation parameter in the problem by a simulated annealing optimization method to obtain a numerical solution of the optimal attenuation parameter that can reduce conservatism; S500, selecting multiple attenuation parameters around the optimal attenuation parameter to obtain different boundaries, generating a new boundary of the reachable set by taking the intersection method, and quantitatively calculating the degree of conservative reduction of the new boundary of the reachable set; Wherein, in the step S100, T p is the basic period n x The state equation of the dimensional time-varying system is expressed as follows: In the formula, Represents n x dimensional state variables; is the identified periodic matrix, υ(t) is the bounded input signal; Wherein, in the step S100, T p The periodic time-varying parameter matrices A(t) and Ψ(t) are expressed as follows: In the formula, A0, A h ,C h n x ×n x Order constant matrix, Ψ0,Ψ h ,Γ h n x ×n u Order constant matrix; Among them, α h (t) and β h (t) are continuous and bounded periodic time-varying functions, and α h (t) = cos(hω p t),β h (t) = sin(hω p t); Wherein, the step S200 includes: S201. Construct piecewise Lyapunov function V m (t) = x T (t)P m (t)x(t)(m=0,1,...,M-1); where P m (t) is a periodic symmetric positive definite matrix and adopts a linear piecewise form. p +mδ,lT p +(m+1)δ) is expressed as: P m (t)=P m +σ(t)P m+1 In the formula, the equal division parameter M represents the number of equal divisions in a cycle, and the width of the sub-interval is σ(t)∈[0,1) is a linear interpolation function; S202, based on the known Lyapunov function reachable set estimation time-varying criterion, the non-convex time-varying function Q m (t)≤0 negative semidefinite problem is transformed into multiple affine time-varying matrix polynomials Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)] negative semidefinite problem; Wherein, in step S202: The known time-varying criterion Q for reaching set estimation based on Lyapunov function m (t) is expressed as follows: In the formula, About V m (t) is the Dini derivative; λ is the decay parameter of the Lyapunov function; Among them, Q m The negative semidefinite form of the quadratic form of (t) is expressed as follows: And the time-varying matrix polynomial Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)] is expressed as follows: In the formula, Λ 1,1,m ,…,Λ 1,2H+1,m and Λ 2,1,m ,…,Λ 2,2H,m They are all constant matrices determined by the system modeling state equation and the Lyapunov function form; Wherein, the step S300 includes: S301, according to the properties of the affine function, a time-varying matrix polynomial is obtained. The necessary and sufficient condition for satisfying the following inequality is that the matrices corresponding to the endpoints of the time-varying parameter value range of the time-varying polynomial are all negative semidefinite; L m [σ(t);α1(t),…,α H (t);β1(t),…,β H (t)]≤0 According to the matrix polynomial Λ m [σ(t);α1(t),…,α H (t); β1(t),…,β H (t)] is continuous and has a time-varying parameter σ(t)∈[0,1), and the matrix polynomial limit set at σ(t)→1 is: S302, according to the lemma of time-varying reachable sets, P m (t) satisfies P m (t) = P m (t+T p )>0, and according to P m (t) Using linear equal interpolation, we can get P m The necessary and sufficient conditions for (t)>0 are: P0>0,P1>0,…,P M-1 >0 S303, through the positive definite matrix P m (t) is scaled, and the following formula is used as the optimization function: In the formula, This is true for m=0,1,...,M-1; According to the constraints Determine the reachable set boundary.

2. The method according to claim 1, characterized in that The step S400 includes: S401, using the cosine sampling function to construct the following attenuation parameters λ0,λ1,…,λ M-1 : Among them, the auxiliary parameter p∈[0,0.5] is set to make all the attenuation parameters λ m >0; S402, obtaining a numerical solution of an optimal attenuation parameter λ for reducing the conservatism of the reachable set by using a simulated annealing optimization algorithm; Among them, the solution process of the simulated annealing algorithm is: Initialize the decay parameter λ = λ (0) After that, the update law of the decay parameter is λ (k+1) =λ (k) +ρ·rand(0,1), where the maximum update step size ρ∈[0,1), rand(0,1) represents the random number generated in the interval [0,1]; when λ=λ (k+1) , substitute into the optimization function to calculate the optimal And calculate Δf= value; Set the current annealing temperature to T, if Δf<0 or Choose to accept λ=λ (k+1) , and set the next iteration temperature to μT, where the decay rate μ∈[0,1); repeat the above steps until the given convergence criterion is reached, and finally obtain the optimal λ value through numerical solution.

3. The method according to claim 2, characterized in that The step S500 includes: S501, select G different λ values ​​near the obtained optimal attenuation parameter λ and substitute them into the optimization problem to obtain G positive definite matrices Constructing new boundaries It is expressed as follows: Among them, the boundary equations L1, L2, ..., L G Respectively by get, is the result of the above different boundary equations L1, L2, ..., L G Find the new boundary obtained by intersection; S502, according to the selected positive definite matrix Obtained by The determined ellipse area S1 or ellipsoid volume V1 is evenly scattered in the space to obtain the area S2 or volume V2 inside the reachable set, and then the efficiency of quantitatively obtaining the exact area or volume of the reachable set is obtained by the intersection method, which is expressed as follows: Where η S It represents the exact efficiency of the reachable set area relative to the original area after using the set intersection method; η V It represents the exact efficiency of the reachable set volume relative to the original volume. 4 . A computer-readable storage medium having program instructions stored thereon, wherein the program instructions, when executed by a processor, implement the method according to claim 1 .

5. A reachable set estimation system for periodic time-varying mechanical systems based on set intersection, characterized in that: include: A computer device comprising a computer-readable storage medium according to claim 4; A periodic time-varying mechanical system, and a system information identification processing module and a plurality of sensors for feeding back system information to the computer device.

Images