CN113110068B - Subspace system identification method and system - Google Patents

Subspace system identification method and system Download PDF

Info

Publication number
CN113110068B
CN113110068B CN202110561327.6A CN202110561327A CN113110068B CN 113110068 B CN113110068 B CN 113110068B CN 202110561327 A CN202110561327 A CN 202110561327A CN 113110068 B CN113110068 B CN 113110068B
Authority
CN
China
Prior art keywords
matrix
list
singular value
value decomposition
output data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110561327.6A
Other languages
Chinese (zh)
Other versions
CN113110068A (en
Inventor
夏元清
高润泽
詹玉峰
翟弟华
戴荔
孙中奇
张金会
闫莉萍
刘坤
郭泽华
崔冰
邹伟东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202110561327.6A priority Critical patent/CN113110068B/en
Publication of CN113110068A publication Critical patent/CN113110068A/en
Application granted granted Critical
Publication of CN113110068B publication Critical patent/CN113110068B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/02—Total factory control, e.g. smart factories, flexible manufacturing systems [FMS] or integrated manufacturing systems [IMS]

Abstract

The invention relates to a subspace system identification method system, which comprises the following steps: initializing four groups of input and output data matrixes, and merging a first input data matrix and a first output data matrix to obtain a first intermediate variable matrix; merging the third input data matrix and the third output data matrix to obtain a second intermediate variable matrix; obtaining a first projection of a row space of the second output data matrix along a row space of the second input data matrix toward a row space of the first intermediate variable matrix; obtaining a second projection of the row space of the fourth output data matrix along the row space of the fourth input data matrix towards the row space of the second intermediate variable matrix; dividing the first projection into column blocks and performing singular value decomposition; aggregating the first singular value decomposition results; performing singular value decomposition on the polymerization result to obtain a second singular value decomposition result; and aggregating the second singular value decomposition result, the first projection and the second projection to obtain a system state space model. The invention improves the processing speed.

Description

Subspace system identification method and system
Technical Field
The present invention relates to the field of system identification technologies, and in particular, to a subspace system identification method and system.
Background
System identification is the process of mathematical modeling of a dynamic system using measured data, and the model established can be used for system analysis, performance monitoring and diagnostics, prediction, optimization, and system design and control. The subspace system identifies the concepts in the aspects of comprehensive system theory, linear algebra, statistics and the like, and directly utilizes the measurement input and output data to identify the system mathematical model. Compared with other system identification methods, the subspace system identification method has the following advantages: the method has the advantages that parameterization is not needed, iterative optimization is not needed, the algorithm only depends on simple and reliable linear algebraic tools such as singular value decomposition, and the like, and the state space model is directly estimated. However, at present, subspace system identification is directly based on a large amount of data, and the calculation speed is limited, so that the subspace system identification is not suitable for processing identification tasks with large scale or real-time constraint conditions. The workflow is an important information processing mode in cloud computing, and a scientific computing task is decomposed into a Directed Acyclic Graph (DAG) structure with a front-back dependency relationship, wherein workflow nodes represent subtasks with small granularity, and edges represent the dependency relationship among the subtasks. And distributing each subtask to distributed computing nodes of the cloud resource pool according to resource requirements, so that the processing efficiency of the scientific computing task is improved. The current subspace system identification research only considers improving the algorithm structure and improving the algorithm performance, and the processing speed is still to be improved.
Disclosure of Invention
The invention aims to provide a subspace system identification method and system, which improve the processing speed.
In order to achieve the purpose, the invention provides the following scheme:
a subspace system identification method includes:
collecting historical input data and historical output data of a controlled system;
according to the historical input data and the historical output data, four groups of input and output data matrixes are initialized, the first group of input and output data matrixes are a first input data Hankel matrix and a first output data Hankel matrix, the second group of input and output data matrixes are a second input data Hankel matrix and a second output data Hankel matrix, the third group of input and output data matrixes are a third input data Hankel matrix and a third output data Hankel matrix, and the fourth group of input and output data matrixes are a fourth input data Hankel matrix and a fourth output data Hankel matrix; the data in the first set of input-output data matrices is early data of the data in the second set of input-output data matrices, and the data in the third set of input-output data matrices is early data of the data in the fourth set of input-output data matrices;
combining the first input data Hankel matrix and the first output data Hankel matrix to obtain a first intermediate variable matrix;
combining the third input data Hankel matrix and the third output data Hankel matrix to obtain a second intermediate variable matrix;
acquiring projection of the line space of the second output data Hankel matrix along the line space of the second input data Hankel matrix to the line space of the first intermediate variable matrix, and defining the projection as a first projection;
acquiring projection of the line space of the fourth output data Hankel matrix along the line space of the fourth input data Hankel matrix to the line space of the second intermediate variable matrix, and defining the projection as a second projection;
dividing the first projection into col matrix blocks by columns; col is an even number;
processing the col matrix blocks in sequence based on a singular value decomposition method to obtain a plurality of first singular value decomposition results;
aggregating every two first singular value decomposition results in the obtained plurality of first singular value decomposition results into one group;
judging whether the number of the obtained aggregation results is 1 or not;
if the number of the obtained aggregation results is 1, outputting the aggregation results;
if the obtained aggregation result is not 1, performing singular value decomposition on the multiple aggregation results to obtain first singular value decomposition results, and returning to the step of aggregating every two first singular value decomposition results in the obtained multiple first singular value decomposition results into a group;
performing singular value decomposition on the polymerization result to obtain a second singular value decomposition result;
and aggregating the second singular value decomposition result, the first projection and the second projection to obtain a system state space model.
Optionally, the aggregating every two first singular value decomposition results in the obtained plurality of first singular value decomposition results into a group specifically includes:
aggregating every two first singular value decomposition results in the obtained multiple first singular value decomposition results into a group by using a BlockMerge function to obtain an aggregation result;
the input of the BlockMerge function is two first singular value decomposition results, and the first decomposition result is U 1 ,∑ 1 ,V 1 And the second decomposition result is U 2 ,∑ 2 ,V 2 The output is the aggregation result U r ,∑ r ,V r ；
Obtaining an aggregation result U through a BlockMerge function r ,∑ r ,V r The process comprises the following steps:
obtain U 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k ；
Obtain U 2 ,∑ 2 ,V 2 Low order approximation result U of 2l ,∑ 2l ,V 2l ；
To U 1k ∑ 1k 、U 2l ∑ 2l Carrying out aggregation and singular value decomposition to obtain U r ,∑ r ,
Figure BDA0003078979620000031
According to the formula
Figure BDA0003078979620000032
To V 1k ,V 2l Carrying out polymerization to obtain V r 。
Optionally, the obtaining U 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k The method specifically comprises the following steps:
obtaining U through DoTruncate 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k ；
Obtaining U through DoTruncate 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k The specific process comprises the following steps:
by the formula k = rank (∑) 1 ) Calculating a singular value order k;
using formula U 1k ＝U 1 (: 1 1 Carrying out low-order approximation processing to obtain U 1k ；
Using the formula Σ 1k ＝∑ 1 (1, k,1 1 Carrying out low-order approximation processing to obtain sigma 1k
Using formula V 1k ＝V 1 (1, k,1 1 Carrying out low-order approximate treatment to obtain V 1k 。
Optionally, the processing, based on the singular value decomposition method, the col matrix blocks in sequence to obtain a plurality of first singular value decomposition results specifically includes:
performing singular value decomposition on the col matrix blocks by using a DoMergeOfBlocks function;
the singular value decomposition of col matrix blocks by using the DoMergeOfBlocks function specifically includes:
preset list l for storing decomposition results U L sigma and l V Wherein, list l U In which the U matrix, list l, of singular value decomposition of each matrix block is stored ∑ In which the sigma matrix of each matrix block subjected to singular value decomposition is stored, list l V A V matrix for performing singular value decomposition on each matrix block is stored in the memory, and the singular value decomposition result of the matrix block is U sigma V T ；
According to the formula Nl = len (l) U ) Obtaining the length Nl of the list;
according to the formula level = ceil (log) 2 Nl), obtaining an iteration layer level;
initializing i =1;
an empty list l is built Ut 、l ∑t And l Vt Use list l U Is 1 of Ut Assign value, list l ∑ Is 1 ∑t Assign value, list l V Is 1 Vt Carrying out assignment;
initializing j =1;
utilizing the BlockMerge function according to a formula U j ,∑ j ,V j ＝BlockMerge(l Ut (j),l ∑t (j),l Vt (j),l Ut (j+1),l ∑t (j+1),l Vt (j + 1)), U is obtained j ,∑ j ,V j (ii) a Wherein l Ut (j) Representation List l Ut The j-th element, l ∑t (j) Representation list l ∑t The j-th element, l Vt (j) Representation list l Vt The jth element.
Will U j Put into List l U In (1), dividing j Put into List l U In (1), V j Put into List l V The preparation method comprises the following steps of (1) performing;
adding 2 to the value of j;
judging whether j is smaller than or equal to the length Nl of the list;
if j is less than or equal to the length Nl of the list, returning to the step of utilizing the BlockMerge function according to the formula U j ,∑ j ,V j ＝BlockMerge(l Ut (j),l ∑t (j),l Vt (j),l Ut (j+1),l ∑t (j+1),l Vt (j + 1)), obtaining U j ,∑ j ,V j ”；
If j is larger than the length Nl of the list, judging whether the length Nl of the list is an odd number;
if the list length Nl is odd, the list l is added Ut Is added to the list l U Will list l ∑t Is added to the list l ∑ Will list l Vt Is added to the list l V ；
If the length Nl of the list is an even number, adding 1 to the value of i;
judging whether i is less than or equal to the iteration layer level;
if i is less than or equal toEqual to the level of the iteration layer, returning to the' built empty list l Ut 、l ∑t And l Vt Use list l U Is 1 Ut Assign value, list l ∑ Is 1 ∑t Assign value, list l V Is 1 of Vt Assigned value'
If i is greater than the iteration layer level, ending the loop;
will list l U 、l ∑ And l V And outputting the result as the first singular value decomposition result.
Optionally, the performing singular value decomposition on the aggregation result to obtain a second singular value decomposition result specifically includes:
according to formula N c = round (n/col' + 0.45) get cut column block number; n represents the number of columns of the aggregation result matrix M to be decomposed, and col' represents the width of a column block after preset splitting;
establishing a list l for storing the cut column blocks A And a list l for storing the decomposition results in each column of blocks U 、l ∑ And l V ；
Cutting matrix M into N columns c Column block, N c Put column blocks into list l A Performing the following steps;
to N c Singular value decomposition is carried out on each column block, and the result is correspondingly stored into a list l U 、l ∑ And l V Performing the following steps;
will l U 、l ∑ And l V The singular value decomposition result of the matrix M is obtained as an input to the domargeofblocks function.
Optionally, the aggregating the second singular value decomposition result, the first projection, and the second projection to obtain a system state space model specifically includes:
cutting the second singular value decomposition results U, S and V, and taking the order n of the singular value S as the system order;
cutting singular value decomposition results U and V according to the system order n to obtain a matrix U 1 And V 1 (ii) a Wherein S, U and V are represented as:
Figure BDA0003078979620000051
U 1 number of columns and V 1 Are all n, U 2 And V 2 For matrices with the remainder discarded, T denotes transpose, S 1 Is a first intermediate matrix, S 2 Is a second intermediate matrix;
computing an expansion matrix Γ i And an expansion matrix Γ i-1 (ii) a Wherein the content of the first and second substances,
Figure BDA0003078979620000052
W 1 representing a weight matrix, Γ i-1 Is gamma i Removing the matrix of the last row;
computing a state estimation sequence X i And X i+1 (ii) a Wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0003078979620000053
representing the pseudo-inverse of the matrix, O i Represents the first projection, O i-1 Representing the second projection to be aggregated;
estimating a sequence X from a state i And X i+1 Obtaining a system state space model; wherein the system state space model is represented by matrices A, B, C and D; matrices A, B, C and D and state estimation sequence X i And X i+1 The relation of (A) is as follows:
Figure BDA0003078979620000054
wherein, Y i Represents X i Corresponding to the output data, U, of the controlled system in the time range i Represents X i Corresponding to the input data vector in the time range.
Optionally, the col value is 10.
The invention also discloses a subspace system identification system, comprising:
the data acquisition module is used for acquiring historical input data and historical output data of a controlled system;
the data initialization module is used for initializing four groups of input and output data matrixes according to the historical input data and the historical output data, wherein the first group of input and output data matrixes are a first input data Hankel matrix and a first output data Hankel matrix, the second group of input and output data matrixes are a second input data Hankel matrix and a second output data Hankel matrix, the third group of input and output data matrixes are a third input data Hankel matrix and a third output data Hankel matrix, and the fourth group of input and output data matrixes are a fourth input data Hankel matrix and a fourth output data Hankel matrix; the data in the first set of input-output data matrices is early data of the data in the second set of input-output data matrices, and the data in the third set of input-output data matrices is early data of the data in the fourth set of input-output data matrices;
the first intermediate variable matrix obtaining module is used for merging the first input data Hankel matrix and the first output data Hankel matrix to obtain a first intermediate variable matrix;
the second intermediate variable matrix obtaining module is used for merging the third input data Hankel matrix and the third output data Hankel matrix to obtain a second intermediate variable matrix;
a first projection obtaining module, configured to obtain a projection of a row space of the second output data Hankel matrix to a row space of the first intermediate variable matrix along the row space of the second input data Hankel matrix, and define the projection as a first projection;
a second projection obtaining module, configured to obtain a projection that a row space of the fourth output data Hankel matrix is projected to a row space of the second intermediate variable matrix along the row space of the fourth input data Hankel matrix, and the projection is defined as a second projection;
the first projection column-by-column division module is used for dividing the first projection column-by-column into col matrix blocks; col is an even number;
the matrix block singular value decomposition module is used for sequentially processing col matrix blocks based on a singular value decomposition method to obtain a plurality of first singular value decomposition results;
the first aggregation module is used for aggregating every two first singular value decomposition results in the obtained plurality of first singular value decomposition results into one group;
the judging module is used for judging whether the number of the obtained aggregation results is 1 or not;
the aggregation result output module is used for outputting the aggregation result if the number of the obtained aggregation results is 1;
the second aggregation module is used for performing singular value decomposition on the multiple aggregation results to obtain first singular value decomposition results if the obtained aggregation result is not 1, and returning to the step of aggregating every two first singular value decomposition results in the multiple obtained first singular value decomposition results into a group;
the aggregation result singular value decomposition module is used for performing singular value decomposition on the aggregation result to obtain a second singular value decomposition result;
and the system state space model obtaining module is used for aggregating the second singular value decomposition result, the first projection and the second projection to obtain a system state space model.
According to the specific embodiment provided by the invention, the invention discloses the following technical effects:
the invention utilizes the historical data of the controlled system to identify the subspace system based on the workflow structure, thereby improving the processing speed of the system identification task.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required in the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art that other drawings can be obtained according to these drawings without creative efforts.
FIG. 1 is a schematic flow chart of a subspace system identification method according to the present invention;
FIG. 2 is a schematic diagram of a subspace system identification system according to the present invention;
FIG. 3 is a schematic diagram of a subspace system identification workflow configuration of the present invention;
FIG. 4 is a schematic view of a club system control configuration of the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in further detail below.
Fig. 1 is a schematic flow chart of a subspace system identification method of the present invention, and as shown in fig. 1, the subspace system identification method includes:
step 101: collecting historical input data and historical output data of a controlled system;
step 102: according to the historical input data and the historical output data, four groups of input and output data matrixes are initialized, the first group of input and output data matrixes are a first input data Hankel matrix and a first output data Hankel matrix, the second group of input and output data matrixes are a second input data Hankel matrix and a second output data Hankel matrix, the third group of input and output data matrixes are a third input data Hankel matrix and a third output data Hankel matrix, and the fourth group of input and output data matrixes are a fourth input data Hankel matrix and a fourth output data Hankel matrix; the data in the first set of input-output data matrices is early data of the data in the second set of input-output data matrices, and the data in the third set of input-output data matrices is early data of the data in the fourth set of input-output data matrices;
step 103: combining the first input data Hankel matrix and the first output data Hankel matrix to obtain a first intermediate variable matrix;
step 104: combining the third input data Hankel matrix and the third output data Hankel matrix to obtain a second intermediate variable matrix;
step 105: acquiring projection of the line space of the second output data Hankel matrix along the line space of the second input data Hankel matrix to the line space of the first intermediate variable matrix, and defining the projection as a first projection;
step 106: acquiring projection of the line space of the fourth output data Hankel matrix along the line space of the fourth input data Hankel matrix to the line space of the second intermediate variable matrix, and defining the projection as a second projection;
step 107: dividing the first projection into col matrix blocks by columns; col is an even number;
step 108: processing the col matrix blocks in sequence based on a singular value decomposition method to obtain a plurality of first singular value decomposition results;
step 109: aggregating every two first singular value decomposition results in the obtained plurality of first singular value decomposition results into one group to obtain an aggregation result;
step 110: judging whether the number of the obtained aggregation results is 1 or not;
step 111: if the obtained aggregation result is not 1, performing singular value decomposition on the multiple aggregation results to obtain a first singular value decomposition result, and returning to the step 109;
step 112: if the number of the obtained aggregation results is 1, outputting the aggregation result
And step 113: performing singular value decomposition on the polymerization result to obtain a second singular value decomposition result;
step 114: and aggregating the second singular value decomposition result, the first projection and the second projection to obtain a system state space model.
And controlling the controlled system according to the system state space model.
In step 101, aggregating every two first singular value decomposition results in the obtained plurality of first singular value decomposition results into a group to obtain an aggregated result, which specifically includes:
aggregating every two first singular value decomposition results in the obtained multiple first singular value decomposition results into a group by using a BlockMerge function to obtain an aggregation result;
the input of the BlockMerge function is two first singular value decomposition results, and the first decomposition result is U 1 ,∑ 1 ,V 1 And the second decomposition result is U 2 ,∑ 2 ,V 2 The output is the aggregation result U r ,∑ r ,V r ；
Obtaining an aggregation result U through a BlockMerge function r ,∑ r ,V r The process comprises the following steps:
obtain U 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k ；
Obtain U 2 ,∑ 2 ,V 2 Low order approximation result U of 2l ,∑ 2l ,V 2l ；
To U 1k ∑ 1k 、U 2l ∑ 2l Carrying out aggregation and singular value decomposition to obtain U r ,∑ r ,
Figure BDA0003078979620000091
According to the formula
Figure BDA0003078979620000092
To V 1k ,V 2l Carrying out polymerization to obtain V r 。
The obtaining U 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k The method specifically comprises the following steps:
obtaining U through DoTruncate 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k ；
Obtaining U through DoTruncate 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k The specific process comprises the following steps:
by the formula k = rank (∑) 1 ) Calculating the order of singular valuek；
Using formula U 1k ＝U 1 (1 1 Carrying out low-order approximation processing to obtain U 1k ；
Using the formula Σ 1k ＝∑ 1 (1, k,1 1 Carrying out low-order approximation processing to obtain sigma 1k
Using formula V 1k ＝V 1 (1 1 Carrying out low-order approximation processing to obtain V 1k 。
Based on a singular value decomposition method, processing col matrix blocks in sequence to obtain a plurality of first singular value decomposition results, specifically comprising:
performing singular value decomposition on col matrix blocks by using a DoMergeOfBlocks function;
the singular value decomposition of col matrix blocks by using the DoMergeOfBlocks function specifically includes:
preset list l for storing decomposition results U 、l ∑ And l V Wherein, list l U In which the U matrix, list l, of singular value decomposition of each matrix block is stored ∑ In which the sigma matrix of each matrix block subjected to singular value decomposition is stored, list l V The V matrix in which singular value decomposition is performed for each matrix block is stored, and the singular value decomposition result of the matrix block is U sigma V T ；
According to the formula Nl = len (l) U ) Obtaining the length Nl of the list;
according to the formula level = ceil (log) 2 Nl), obtaining an iteration layer level;
initializing i =1;
an empty list l is built Ut 、l ∑t And l Vt Use list l U Is 1 of Ut Assign value, list l ∑ Is 1 ∑t Assign value, list l V Is 1 Vt Carrying out assignment;
initializing j =1;
using the BlockMerge function according to a formula U j ,∑ j ,V j ＝BlockMerge(l Ut (j),l ∑t (j),l Vt (j),l Ut (j+1),l ∑t (j+1),l Vt (j + 1)), obtaining U j ,∑ j ,V j (ii) a Wherein l Ut (j) Representation list l Ut The j-th element, l ∑t (j) Representation list l ∑t The j-th element, l Vt (j) Representation list l Vt The j-th element, l Ut (j + 1) denotes a list l Ut J +1 th element, l ∑t (j + 1) denotes a list l ∑t Element j +1, l Vt (j + 1) indicates a list l Vt The j +1 th element.
Will U j Put into List l U In, will ∑ j Put into List l U In (1), V j Put into List l V Performing the following steps;
adding 2 to the value of j;
judging whether j is smaller than or equal to the length Nl of the list;
if j is less than or equal to the length Nl of the list, returning to the step of utilizing the BlockMerge function according to the formula U j ,∑ j ,V j ＝BlockMerge(l Ut (j),l ∑t (j),l Vt (j),l Ut (j+1),l ∑t (j+1),l Vt (j + 1)), obtaining U j ,∑ j ,V j ”；
If j is larger than the length Nl of the list, judging whether the length Nl of the list is an odd number;
if the length Nl of the list is odd, the list l is divided into two parts Ut Is added to the list l U Will list l ∑t Is added to the list l ∑ Will list l Vt Is added to the list l V ；
If the list length Nl is an even number, adding 1 to the value of i;
judging whether i is less than or equal to the iteration layer level;
if i is less than or equal to the level of the iteration layer, returning to the' built empty list l Ut 、l ∑t And l Vt Use list l U Is 1 of Ut Assign value, list l ∑ Is 1 of ∑t Assign value, list l V Is 1 Vt Assigned value'
If i is greater than the iteration layer level, ending the loop;
will list l U 、l ∑ And l V And outputting as a first singular value decomposition result.
The singular value decomposition of the aggregation result to obtain a second singular value decomposition result specifically includes:
according to formula N c = round (n/col' + 0.45) obtain cut column block number; n represents the column number of the aggregation result matrix M to be decomposed, and col' represents the width of a column block after preset splitting;
building a list l for storing the cut column blocks A And a list l for storing the decomposition results in each column of blocks U 、l ∑ And l V ；
Cutting matrix M into N columns c Column block, N c Put column blocks into list l A Performing the following steps;
to N c Singular value decomposition is carried out on each column block, and the results are correspondingly stored in a list l U 、l ∑ And l V The preparation method comprises the following steps of (1) performing;
will l U 、l ∑ And l V The singular value decomposition result of the matrix M is obtained as an input to the domargeofblocks function.
The aggregating the second singular value decomposition result, the first projection and the second projection to obtain a system state space model specifically includes:
cutting the second singular value decomposition result U, S and V, and taking the order n of the singular value S as a system order;
cutting singular value decomposition results U and V according to the system order n to obtain a matrix U 1 And V 1 (ii) a Wherein S, U and V are represented as:
Figure BDA0003078979620000111
U 1 number of columns and V 1 Are all n, U 2 And V 2 For matrices with the remainder discarded, T denotes transpose, S 1 Is a first intermediate matrix, S 2 Is a second intermediate matrix;
calculating the expansion matrix gamma i And an expansion matrix Γ i-1 (ii) a Wherein the content of the first and second substances,
Figure BDA0003078979620000112
Γ i-1 is gamma i Matrix of last row is removed, W 1 Representing a weight matrix, W 1 The superscript of-1 represents the inversion;
computing a state estimation sequence X i And X i+1 (ii) a Wherein the content of the first and second substances,
Figure BDA0003078979620000113
representing the pseudo-inverse of the matrix, O i Representing said first projection, O i-1 Representing the second projection to be aggregated;
estimating a sequence X from a state i And X i+1 Obtaining a system state space model; wherein the system state space model is represented by matrices A, B, C and D; matrices A, B, C and D and state estimation sequence X i And X i+1 The relation of (A) is as follows:
Figure BDA0003078979620000114
wherein Y is i Represents X i Corresponding to the output data, U, of the controlled system in the time range i Represents X i Corresponding to the input data vector in the time range. Y is i ＝[y(N),y(N+1),...,y(N+j-1)] T ，U i ＝[u(N),u(N+1),...,x(N+j-1)] T
ABCD represents a state space model matrix. The state space model matrix is:
x(k+1)＝Ax(k)+Bu(k)；
y(k)＝Cx(k)+Du(k)；
where k is the time, u (k) is the input variable, x (k) is the state variable (hidden in the system, not directly observed, visible as the intermediate variable), and y (k) is the output variable.
Through the iteration of the state space model matrix, the system state can be continuously updated, and the system model is obtained.
The system identification is to calculate the four ABCD matrixes, and then the system model can be obtained.
To further illustrate the application of the present invention, a subspace system identification method is described as an embodiment of subspace system identification for a club system. As shown in fig. 4, the club system comprises a control computer, a servo driver, a club mechanical body (mechanical part) and a sensor, which form a closed loop system, and the control is to make the ball reach a designated position on the cross bar and stop and keep the ball still. At a certain sampling time t, the input variable of the club system is u (t), representing the gear angle theta, and the detected output variable is
Figure BDA0003078979620000121
Wherein y is 1 (t),y 2 (t) represents the position and velocity of the pellet, respectively. The method comprises the steps of identifying and obtaining a system state space model between input and output, analyzing and designing a control algorithm according to an identified mathematical model to realize more accurate control on a controlled system, setting parameters N, j, enabling a club system to run at 2N + j-1 sampling moments, collecting input data and output data of the club system in the period of time, and further processing the following steps:
establishing an input and output data Hankel matrix:
first input data Hankel matrix:
Figure BDA0003078979620000122
second input data Hankel matrix:
Figure BDA0003078979620000123
first output data Hankel matrix:
Figure BDA0003078979620000124
second output data Hankel matrix:
Figure BDA0003078979620000125
third input data Hankel matrix:
Figure BDA0003078979620000126
fourth input data Hankel matrix:
Figure BDA0003078979620000131
third output data Hankel matrix:
Figure BDA0003078979620000132
fourth output data Hankel matrix:
Figure BDA0003078979620000133
then, performing subspace system identification based on the workflow structure:
fig. 3 is a subspace system identification workflow structure, i.e. identification work subtask sequence, for the case of col =10, where col means parallelism. Table 1 describes the subspace system identification workflow task functions and the inter-task dependencies in this example.
TABLE 1 subspace System identification workflow task function and dependency description
Mirror image Serial number Hierarchy level Function(s)
Task ini - - Initializing, creating historical data and preparing Hankel matrix
Task A Task 1 1 Generating O for singular value decomposition tasks i
Task B Task 2 1 Generating O for an expanded matrix computation task i-1
Task C Task 3 -Task 12 2 Is O i Each matrix block of (a) performs a base singular value decomposition operation
Task D Task 13 -Task 19 3，4 Aggregating the results of the parent task; performing a new base singular value decomposition
Task E Task 20 5 Aggregating the results of all tasks; solving a system state space model
In this example, there are 5 task image types in the complete workflow, except for the initial preparation task. For example, mirror image Task C Can be multiplexed 10 times with different input data in the second layer of the workflow. The task functions and relationships are described as follows:
(1) Initial Task ini ，Task ini Is responsible for initializing, creating historical data, preparing the Hankel matrix for the workflow (i.e., creating U in this embodiment) p ,Y p ,U f ,Y f ,U pp ,Y pp ,U fm ,Y fm And is finished).
(2) Based on Hankel matrix [ U ] p ,Y p ,U f ,Y f ]，Task 1 Executing W p ＝[U p ,Y p ]And calculating the oblique projection O i ；
Figure BDA0003078979620000141
Calculated result O i Is divided into 10 column blocks by column, and the 10 column blocks are transmitted to the Task in sequence 3 -Task 12 All of O is i Is passed to Task 20 。
(3) Similarly, based on Hankel matrix [ U ] pp ,Y pp ,U fm ,Y fm ]，Task 2 Calculating the oblique projection O i-1 ；
Figure BDA0003078979620000142
Then adding O i-1 To Task 20 ；
(4)Task 1 The generated 10 column blocks are respectively transmitted to the Task 3 -Task 12 . In each Task, the singular value decomposition is performed on each column block using the function domargeofblocks, and the decomposition results are passed to the Task in the relationship shown in fig. 3 13 -Task 19 (ii) a Here the singular value decomposition is a substantially linear algebraic tool.
(5)Task 13 -Task 19 Using the function blockMerge to the parent Task result (Task) 3 -Task 12 ) Performing aggregation, performing singular value decomposition again on the aggregated result, and transferring the decomposition result to Task according to the relation shown in FIG. 3 20 (ii) a Here the singular value decomposition is a substantially linear algebraic tool.
(6) Finally, the egress Task 20 Collecting all parent tasks (Task) 1-2 ，Task 17-19 ) And aggregating the results, and solving to obtain a state space model of the identified system, which specifically comprises the following steps:
s1, merging Task based on function parallelSVdbyCols 17 -Task 19 To obtain O i The singular value decomposition results U, S, V.
And S2, cutting a singular value decomposition result, and taking the order n of the singular value S as a system order.
S3, cutting the U and V obtained by singular value decomposition according to order n to obtain a matrix U 1 、V 1 。
Figure BDA0003078979620000143
Wherein, U 1 Number of rows, V 1 Number of lines n, U 2 、V 2 The remaining discarded matrices.
S4, calculating an expansion matrix gamma i And gamma i-1 ：
Figure BDA0003078979620000151
Γ i-1 ＝ i Γ
Wherein, the first and the second end of the pipe are connected with each other, i Γrepresents gamma i The matrix of the last row is removed.
S5, according to Task 1-2 Output result of (1) O i And O i-1 Calculating a state estimation sequence X i And X i+1 ；
Figure BDA0003078979620000152
Figure BDA0003078979620000153
Wherein
Figure BDA0003078979620000154
Representing the pseudo-inverse of the matrix.
S6, solving a linear equation to obtain system state space matrixes A, B, C and D;
Figure BDA0003078979620000155
note: the above S4 and S5 correspond to a partial operation of the function ParallelSVdbyCols. Because the workflow structure is optimized, all the circulation in the function ParallelSVdbyCols is not finished, but the residual flow of the S1 singular value decomposition stage is merged to the export Task 20 In, i.e. as in FIG. 2 Task 17 、Task 18 And Task 19 After singular value decomposition, two aggregation are performed, and finally an aggregation result is output, for example: task 18 And Task 19 After singular value decomposition, the singular value decomposition is carried out and the singular value decomposition is aggregated into a result, and the result and the Task are combined 17 After the singular value decomposition, the singular value decomposition is respectively carried out and then the singular value decomposition is carried out, and a polymerization result is obtained, wherein the process is carried out in Task 20 Is carried out in (1).
The process of establishing the four basis functions of the identification method is described below.
Designing four functions of parallelSVDbyCols, doMergeOfBlocks, blockMerge and DoTruncate. In S3, the singular value decomposition stage in S1 is decomposed into a workflow structure based on the four functions, and the workflow structure is finally used as a part of the identification workflow of the complete subspace system.
Singular value decomposition: a = USV T . A is the decomposed matrix, T is superscripted as transposed, and U, S (Σ), and V are the three decomposed matrices.
1. Designing a parallelSVdbyCols function:
inputting: after the M × n matrix M is decomposed, the column block width col is split.
And (3) outputting: singular value decomposition results
Figure BDA0003078979620000156
a1, calculating the number of cutting column blocks: n is a radical of c ＝round(n/col+0.45)。
a2, establishing a list for storing column blocks and singular value decomposition results of the column blocks, and dividing l A And (6) filling.
l A ＝list()；l U ＝list()；l ∑ ＝list()；l V ＝list()；
Where list () is a programming language notation, meaning that an empty list is created, l A Is a list for storing the cut column blocks, l U ,l ∑ ,l V The blocks in each column store a list of decomposition results.
and a3, performing singular value decomposition on each column of blocks, and storing the obtained result into an a2 list.
Note: here singular value decomposition into substantially linear algebraic tools
a4, calling a DoMergeOfBlocks function, and returning an M singular value decomposition result:
Figure BDA0003078979620000161
2. designing a DoMergeOfBlocks function:
inputting: list l of singular value decomposition results for each column of blocks U ,l ∑ ,l V 。
And (3) outputting: singular value decomposition results
Figure BDA0003078979620000162
a1, calculating the length Nl of the list: nl = len (l) U )。
a2, calculating iteration layer level: level = ceil (log) 2 Nl)。
and a3, establishing a double-layer cycle, and solving the M singular value decomposition result.
a3.1, initialization i =1.
a3.2, establishing a new empty list l Ut ,l ∑t ,l Vt Using the input list l U ,l ∑ ,l V And setting the input list to be null for the corresponding assignment:
l Ut ＝l U ,l ∑t ＝l ∑ ,l Vt ＝l V ；
l U ＝list()；l ∑ ＝list()；l V ＝list()。
a3.3, initialization j =1.
a3.4, calling BlockMerge, returning an input two-item data aggregation result, and adding the result into the list l Ut ,l ∑t ,l Vt The method comprises the following steps:
U j ,∑ j ,V j ＝BlockMerge(l Ut (j),l ∑t (j),l Vt (j),l Ut (j+1),l ∑t (j+1),l Vt (j+1))；
l U +＝U j ；l ∑ +＝∑ j ；l V +＝V j 。
a3.5, set j = j +2.
a3.6, judging that j is less than or equal to Nl, if so, returning to the step a3.4, and if not, entering the step a3.7.
and a3.7, judging whether Nl is odd, if so, entering a step a3.8, and if not, entering a step a3.9.
a3.8, will list l Ut ,l ∑t ,l Vt The last entry is added to the list l U ,l ∑ ,l V In (1).
a3.9, setting i = i +1.
and a3.10, judging that i is less than or equal to level, if so, returning to the step a3.2, and if not, ending the circulation.
a4, returning the M singular value decomposition result
Figure BDA0003078979620000171
3. Designing a BlockMerge function:
inputting: two groups of singular value decomposition results U 1 ,∑ 1 ,V 1 ,U 2 ,∑ 2 ,V 2 。
And (3) outputting: results after polymerization U r ,∑ r ,V r 。
a1, calling DoTruncate and returning U 1 ,∑ 1 ,V 1 The low order approximation:
U 1k ,∑ 1k ,V 1k ＝DoTruncate(U 1 ,∑ 1 ,V 1 )。
a2, calling DoTruncate and returning U 2 ,∑ 2 ,V 2 The low order approximation:
U 2l ,∑ 2l ,V 2l ＝DoTruncate(U 2 ,∑ 2 ,V 2 )。
a3, to U 1k ∑ 1k 、U 2l ∑ 2l Polymerization and singular value decomposition:
Figure BDA0003078979620000172
a4, pair V 1k ,V 2l Carrying out polymerization:
Figure BDA0003078979620000173
a5, returning a result U after aggregation r ,∑ r ,V r 。
4. Designing a DoTruncate function:
inputting: and decomposing the results U, sigma and V of the singular value matrix.
And (3) outputting: low order approximation result U k ,∑ k ,V k 。
a1, calculating a singular value order k: k = rank (∑).
a2, respectively carrying out low-order approximation processing on U, sigma and V:
U k ＝U(:,1:k)；∑ k ＝∑(1:k,1:k)；V k ＝V(1:k,1:k)。
a3, returning a low-order approximate resultU k ,∑ k ,V k 。
Fig. 2 is a schematic structural diagram of a subspace system identification system of the present invention, and as shown in fig. 2, the subspace system identification system includes:
the data acquisition module 201 is used for acquiring historical input data and historical output data of a controlled system;
the data initialization module 202 is configured to initialize four groups of input and output data matrices according to the historical input data and the historical output data, where the first group of input and output data matrices are a first input data Hankel matrix and a first output data Hankel matrix, the second group of input and output data matrices are a second input data Hankel matrix and a second output data Hankel matrix, the third group of input and output data matrices are a third input data Hankel matrix and a third output data Hankel matrix, and the fourth group of input and output data matrices are a fourth input data Hankel matrix and a fourth output data Hankel matrix; the data in the first set of input-output data matrices is early data of the data in the second set of input-output data matrices, and the data in the third set of input-output data matrices is early data of the data in the fourth set of input-output data matrices;
a first intermediate variable matrix obtaining module 203, configured to combine the first input data Hankel matrix and the first output data Hankel matrix to obtain a first intermediate variable matrix;
a second intermediate variable matrix obtaining module 204, configured to combine the Hankel matrix of the third input data and the Hankel matrix of the third output data to obtain a second intermediate variable matrix;
a first projection obtaining module 205, configured to obtain a projection that a row space of the second output data Hankel matrix is projected to a row space of the first intermediate variable matrix along the row space of the second input data Hankel matrix, and the projection is defined as a first projection;
a second projection obtaining module 206, configured to obtain a projection that the row space of the fourth output data Hankel matrix is projected to the row space of the second intermediate variable matrix along the row space of the fourth input data Hankel matrix, and define the projection as a second projection;
a first projection column-wise division module 207, configured to divide the first projection column-wise into col matrix blocks; col is an even number;
a matrix block singular value decomposition module 208, configured to sequentially process col matrix blocks based on a singular value decomposition method, so as to obtain a plurality of first singular value decomposition results;
a first aggregation module 209, configured to aggregate every two first singular value decomposition results in the obtained plurality of first singular value decomposition results into a group;
a judging module 210, configured to judge whether the obtained aggregation result number is 1;
an aggregation result output module 212, configured to output an aggregation result if the obtained number of aggregation results is 1;
the second aggregation module 211, if the obtained aggregation result is not 1, is configured to perform singular value decomposition on the multiple aggregation results to obtain first singular value decomposition results, and return to the step "to aggregate every two first singular value decomposition results in the obtained multiple first singular value decomposition results into one group";
an aggregation result singular value decomposition module 212, configured to perform singular value decomposition on the aggregation result to obtain a second singular value decomposition result;
a system state space model obtaining module 213, configured to aggregate the second singular value decomposition result, the first projection, and the second projection, and obtain a system state space model.
According to the specific embodiment provided by the invention, the invention discloses the following technical effects:
the invention utilizes the historical data of the controlled system to carry out the subspace system identification based on the workflow structure, thereby improving the processing speed of the system identification task.
The embodiments in the present description are described in a progressive manner, each embodiment focuses on differences from other embodiments, and the same and similar parts among the embodiments are referred to each other. For the system disclosed by the embodiment, the description is relatively simple because the system corresponds to the method disclosed by the embodiment, and the relevant points can be referred to the method part for description.
The principles and embodiments of the present invention have been described herein using specific examples, which are provided only to help understand the method and the core concept of the present invention; meanwhile, for a person skilled in the art, according to the idea of the present invention, the specific embodiments and the application range may be changed. In view of the above, the present disclosure should not be construed as limiting the invention.

Claims (5)

1. A subspace system identification method, comprising:
collecting historical input data and historical output data of a controlled system;
according to the historical input data and the historical output data, four groups of input and output data matrixes are initialized, the first group of input and output data matrixes are a first input data Hankel matrix and a first output data Hankel matrix, the second group of input and output data matrixes are a second input data Hankel matrix and a second output data Hankel matrix, the third group of input and output data matrixes are a third input data Hankel matrix and a third output data Hankel matrix, and the fourth group of input and output data matrixes are a fourth input data Hankel matrix and a fourth output data Hankel matrix; the data in the first set of input-output data matrices is early data of the data in the second set of input-output data matrices, and the data in the third set of input-output data matrices is early data of the data in the fourth set of input-output data matrices;
combining the first input data Hankel matrix and the first output data Hankel matrix to obtain a first intermediate variable matrix;
combining the third input data Hankel matrix and the third output data Hankel matrix to obtain a second intermediate variable matrix;
acquiring a projection of the line space of the second output data Hankel matrix along the line space of the second input data Hankel matrix to the line space of the first intermediate variable matrix, and defining the projection as a first projection;
acquiring projection of the line space of the fourth output data Hankel matrix along the line space of the fourth input data Hankel matrix to the line space of the second intermediate variable matrix, and defining the projection as a second projection;
dividing the first projection into col matrix blocks by columns; col is an even number;
processing col matrix blocks in sequence based on a singular value decomposition method to obtain a plurality of first singular value decomposition results;
aggregating every two first singular value decomposition results in the obtained plurality of first singular value decomposition results into a group;
judging whether the number of the obtained aggregation results is 1 or not;
if the obtained number of the aggregation results is 1, outputting the aggregation results;
if the obtained aggregation result is not 1, performing singular value decomposition on the multiple aggregation results to obtain first singular value decomposition results, and returning to the step of aggregating every two first singular value decomposition results in the obtained multiple first singular value decomposition results into a group;
performing singular value decomposition on the polymerization result to obtain a second singular value decomposition result;
aggregating the second singular value decomposition result, the first projection and the second projection to obtain a system state space model;
the ball arm system includes control computer, servo driver, ball arm mechanical body and sensor, and the control computer, servo driver, ball arm mechanical body and sensor form a closed loop system, and the control is aimed at making small ball reach a certain defined position on the cross bar, and stopping and holding it, and at a certain sampling time t, the input variable of ball arm system is u (t), and represents gear angle theta, and the detected output variable is u (t)
Figure FDA0003942959790000021
Wherein y is 1 (t),y 2 (t) represents the position and velocity of the pellet, respectively;
enabling the club system to operate 2N + j-1 sampling moments, and collecting input data and output data of the club system in the 2N + j-1 sampling moments;
establishing an input and output data Hankel matrix:
first input data Hankel matrix:
Figure FDA0003942959790000022
second input data Hankel matrix:
Figure FDA0003942959790000023
first output data Hankel matrix:
Figure FDA0003942959790000024
second output data Hankel matrix:
Figure FDA0003942959790000025
third input data Hankel matrix:
Figure FDA0003942959790000031
fourth input data Hankel matrix:
Figure FDA0003942959790000032
third output data Hankel matrix:
Figure FDA0003942959790000033
fourth output data Hankel matrix:
Figure FDA0003942959790000034
performing subspace system identification based on a workflow structure;
aggregating every two singular value decomposition results in the obtained plurality of singular value decomposition results into a group specifically includes:
aggregating every two first singular value decomposition results in the obtained multiple first singular value decomposition results into a group by using a BlockMerge function to obtain an aggregation result;
the input of the BlockMerge function is two first singular value decomposition results, and the first decomposition result is U 1 ,∑ 1 ,V 1 And the second decomposition result is U 2 ,∑ 2 ,V 2 The output is the aggregation result U r ,∑ r ,V r ；
Obtaining an aggregation result U through a BlockMerge function r ,∑ r ,V r The process comprises the following steps:
obtain U 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k ；
Obtain U 2 ,∑ 2 ,V 2 Low order approximation result U of 2l ,∑ 2l ,V 2l ；
To U 1k ∑ 1k 、U 2l ∑ 2l Performing polymerization and singular value decomposition to obtain
Figure FDA0003942959790000035
According to the formula
Figure FDA0003942959790000036
To V 1k ,V 2l Polymerization is carried out to obtain V r ；
The method for processing col matrix blocks in sequence based on singular value decomposition obtains a plurality of first singular value decomposition results, and specifically comprises the following steps:
performing singular value decomposition on col matrix blocks by using a DoMergeOfBlocks function;
the singular value decomposition of col matrix blocks by using the DoMergeOfBlocks function specifically includes:
presetting a list l for storing decomposition results U 、l ∑ And l V Wherein, list l U A list l of U matrices in which the singular value decomposition is performed for each matrix block is stored ∑ Storing the sigma matrix of each matrix block for singular value decomposition, and listing V The V matrix in which singular value decomposition is performed for each matrix block is stored, and the singular value decomposition result of the matrix block is U sigma V T ；
According to the formula Nl = len (l) U ) Obtaining the length Nl of the list;
according to the formula level = ceil (log) 2 Nl), obtaining an iteration layer level;
initializing i =1;
an empty list l is built Ut 、l ∑t And l Vt Use list l U Is 1 of Ut Assign value, list l ∑ Is 1 ∑t Assign value, list l V Is 1 of Vt Carrying out assignment;
initializing j =1;
using the BlockMerge function according to a formula U j ,∑ j ,V j ＝BlockMerge(l Ut (j),l ∑t (j),l Vt (j),l Ut (j+1),l ∑t (j+1),l Vt (j + 1)), obtaining U j ,∑ j ,V j (ii) a Wherein l Ut (j) Representation list l Ut The j-th element, l ∑t (j) Representation list l ∑t The j-th element, l Vt (j) Representation list l Vt The jth element;
will U j Put into List l U In, will ∑ j Put into List l U In (1), V j Put into List l V Performing the following steps;
adding 2 to the value of j;
judging whether j is smaller than or equal to the length Nl of the list;
if j is less than or equal to the length Nl of the list, then returning to the step of utilizing the BlockMerge function according to a formula U j ,∑ j ,V j ＝BlockMerge(l Ut (j),l ∑t (j),l Vt (j),l Ut (j+1),l ∑t (j+1),l Vt (j + 1)), obtaining U j ,∑ j ,V j ”；
If j is larger than the length Nl of the list, judging whether the length Nl of the list is an odd number;
if the list length Nl is odd, the list l is added Ut Is added to the list l U Will list l ∑t Is added to the list l ∑ Will list l Vt Is added to the list l V ；
If the length Nl of the list is an even number, adding 1 to the value of i;
judging whether i is less than or equal to the iteration layer level;
if i is less than or equal to the level of the iteration layer, returning to the' built empty list l Ut 、l ∑t And l Vt Use list l U Is 1 Ut Assign value, list l ∑ Is 1 of ∑t Assign value, list l V Is 1 of Vt Assigning values ";
if i is greater than the iteration layer level, ending the loop;
will list l U 、l ∑ And l V Output as a first singular value decomposition result;
the singular value decomposition of the aggregation result to obtain a second singular value decomposition result specifically includes:
according to formula N c = round (n/col' + 0.45) obtain cut column block number; n represents the column number of the aggregation result matrix M to be decomposed, and col' represents the width of a column block after preset splitting;
establishing a list l for storing the cut column blocks A And a list l for storing the decomposition results in each column of blocks U 、l ∑ And l V ；
Cutting matrix M into N columns c Column block, N c Put a column block into the list l A Performing the following steps;
to N c Singular value decomposition is carried out on each column block, and the result is correspondingly stored into a list l U 、l ∑ And l V Performing the following steps;
will l U 、l ∑ And l V The singular value decomposition result of the matrix M is obtained as an input to the domargeofblocks function.
2. The subspace system identification method of claim 1, wherein the obtaining U 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k The method specifically comprises the following steps:
obtaining U through DoTruncate 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k ；
Obtaining U through DoTruncate 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k The specific process comprises the following steps:
by the formula k = rank (∑) 1 ) Calculating a singular value order k;
using formula U 1k ＝U 1 (: 1 1 Carrying out low-order approximation processing to obtain U 1k ；
Using the formula ∑ 1k ＝∑ 1 (1 1 Carrying out low-order approximation processing to obtain sigma 1k ；
Using formula V 1k ＝V 1 (1, k,1 1 Carrying out low-order approximation processing to obtain V 1k 。
3. The subspace system identification method according to claim 1, wherein the aggregating the second singular value decomposition result, the first projection, and the second projection to obtain a system state space model specifically includes:
cutting the second singular value decomposition results U, S and V, and taking the order n of the singular value S as the system order;
cutting singular value decomposition results U and V according to the system order n to obtain a matrix U 1 And V 1 (ii) a Wherein S, U and V are represented as:
Figure FDA0003942959790000051
U 1 number of columns and V 1 Are all n, U 2 And V 2 Is discarded for the remainderT denotes transpose, S 1 Is a first intermediate matrix, S 2 Is a second intermediate matrix;
computing an expansion matrix Γ i And an expansion matrix Γ i-1 (ii) a Wherein the content of the first and second substances,
Figure FDA0003942959790000061
W 1 representing a weight matrix, Γ i-1 Is gamma i Removing the matrix of the last row;
computing a state estimation sequence X i And X i+1 (ii) a Wherein the content of the first and second substances,
Figure FDA0003942959790000062
Figure FDA0003942959790000063
representing the pseudo-inverse of the matrix, O i Representing said first projection, O i-1 Representing the second projection as being aggregated;
estimating a sequence X from a state i And X i+1 Obtaining a system state space model; wherein the system state space model is represented by matrices A, B, C and D; matrices A, B, C and D and state estimation sequence X i And X i+1 The relation of (A) is as follows:
Figure FDA0003942959790000064
wherein Y is i Represents X i Corresponding to the output data, U, of the controlled system in the time range i Represents X i Corresponding to the input data vector in the time range.
4. The subspace system identification method of claim 1, wherein the value of col is 10.
5. A subspace system identification system, comprising:
the data acquisition module is used for acquiring historical input data and historical output data of a controlled system;
the data initialization module is used for initializing four groups of input and output data matrixes according to the historical input data and the historical output data, wherein the first group of input and output data matrixes are a first input data Hankel matrix and a first output data Hankel matrix, the second group of input and output data matrixes are a second input data Hankel matrix and a second output data Hankel matrix, the third group of input and output data matrixes are a third input data Hankel matrix and a third output data Hankel matrix, and the fourth group of input and output data matrixes are a fourth input data Hankel matrix and a fourth output data Hankel matrix; the data in the first set of input-output data matrices is early data of the data in the second set of input-output data matrices, and the data in the third set of input-output data matrices is early data of the data in the fourth set of input-output data matrices;
the first intermediate variable matrix obtaining module is used for merging the first input data Hankel matrix and the first output data Hankel matrix to obtain a first intermediate variable matrix;
the second intermediate variable matrix obtaining module is used for merging the third input data Hankel matrix and the third output data Hankel matrix to obtain a second intermediate variable matrix;
a first projection obtaining module, configured to obtain a projection of a row space of the second output data Hankel matrix to a row space of the first intermediate variable matrix along the row space of the second input data Hankel matrix, and define the projection as a first projection;
a second projection obtaining module, configured to obtain a projection that a row space of the fourth output data Hankel matrix is projected to a row space of the second intermediate variable matrix along the row space of the fourth input data Hankel matrix, and the projection is defined as a second projection;
the first projection column-by-column division module is used for dividing the first projection column-by-column into col matrix blocks; col is an even number;
the matrix block singular value decomposition module is used for sequentially processing col matrix blocks based on a singular value decomposition method to obtain a plurality of first singular value decomposition results;
the first aggregation module is used for aggregating every two first singular value decomposition results in the obtained plurality of first singular value decomposition results into one group;
the judging module is used for judging whether the number of the obtained aggregation results is 1 or not;
the aggregation result output module is used for outputting the aggregation result if the number of the obtained aggregation results is 1;
the second aggregation module is used for performing singular value decomposition on the multiple aggregation results to obtain first singular value decomposition results if the obtained aggregation result is not 1, and returning to the step of aggregating every two first singular value decomposition results in the multiple obtained first singular value decomposition results into a group;
the aggregation result singular value decomposition module is used for performing singular value decomposition on the aggregation result to obtain a second singular value decomposition result;
a system state space model obtaining module, configured to aggregate the second singular value decomposition result, the first projection, and the second projection to obtain a system state space model;
the ball arm system includes control computer, servo driver, ball arm mechanical body and sensor, and the control computer, servo driver, ball arm mechanical body and sensor form a closed loop system, and the control is aimed at making small ball reach a certain defined position on the cross bar, and stopping and holding it, and at a certain sampling time t, the input variable of ball arm system is u (t), and represents gear angle theta, and the detected output variable is u (t)
Figure FDA0003942959790000081
Wherein y is 1 (t),y 2 (t) represents the position and velocity of the pellet, respectively;
enabling the club system to run at 2N + j-1 sampling moments, and collecting input data and output data of the club system at the 2N + j-1 sampling moments;
establishing an input and output data Hankel matrix:
first input data Hankel matrix:
Figure FDA0003942959790000082
second input data Hankel matrix:
Figure FDA0003942959790000083
first output data Hankel matrix:
Figure FDA0003942959790000084
second output data Hankel matrix:
Figure FDA0003942959790000085
third input data Hankel matrix:
Figure FDA0003942959790000086
fourth input data Hankel matrix:
Figure FDA0003942959790000091
third output data Hankel matrix:
Figure FDA0003942959790000092
fourth output data Hankel matrix:
Figure FDA0003942959790000093
performing subspace system identification based on a workflow structure;
aggregating every two first singular value decomposition results in the obtained plurality of first singular value decomposition results into a group, specifically including:
aggregating every two first singular value decomposition results in the obtained multiple first singular value decomposition results into a group by using a BlockMerge function to obtain an aggregation result;
the input of the BlockMerge function is two first singular value decomposition results, and the first decomposition result is U 1 ,∑ 1 ,V 1 And the second decomposition result is U 2 ,∑ 2 ,V 2 The output is the aggregation result U r ,∑ r ,V r ；
Obtaining an aggregation result U through a BlockMerge function r ,∑ r ,V r The process comprises the following steps:
obtain U 1 ,∑ 1 ,V 1 Low order approximation result U of 1k ,∑ 1k ,V 1k ；
Obtain U 2 ,∑ 2 ,V 2 Low order approximation result U of 2l ,∑ 2l ,V 2l ；
To U 1k ∑ 1k 、U 2l ∑ 2l Performing polymerization and singular value decomposition to obtain
Figure FDA0003942959790000094
According to the formula
Figure FDA0003942959790000095
To V 1k ,V 2l Carrying out polymerization to obtain V r ；
The method for processing col matrix blocks in sequence based on singular value decomposition obtains a plurality of first singular value decomposition results, and specifically comprises the following steps:
performing singular value decomposition on col matrix blocks by using a DoMergeOfBlocks function;
the singular value decomposition of col matrix blocks by using the DoMergeOfBlocks function specifically includes:
preset list l for storing decomposition results U 、l ∑ And l V Wherein, list l U In which the U matrix, list l, of singular value decomposition of each matrix block is stored ∑ In which the sigma matrix of each matrix block subjected to singular value decomposition is stored, list l V In which each matrix block is stored for singular value decompositionThe singular value decomposition result of the matrix block is U sigma V T ；
According to the formula Nl = len (l) U ) Obtaining the length Nl of the list;
according to the formula level = ceil (log) 2 Nl), obtaining an iteration layer level;
initializing i =1;
an empty list l is built Ut 、l ∑t And l Vt Using the list l U Is 1 of Ut Assign value, list l ∑ Is 1 ∑t Assign value, list l V Is 1 Vt Carrying out assignment;
initializing j =1;
utilizing the BlockMerge function according to a formula U j ,∑ j ,V j ＝BlockMerge(l Ut (j),l ∑t (j),l Vt (j),l Ut (j+1),l ∑t (j+1),l Vt (j + 1)), obtaining U j ,∑ j ,V j (ii) a Wherein l Ut (j) Representation list l Ut The j-th element, l ∑t (j) Representation list l ∑t The j-th element, l Vt (j) Representation list l Vt The jth element;
will U j Put into List l U In, will ∑ j Put into List l U In (1), V j Put into List l V The preparation method comprises the following steps of (1) performing;
adding 2 to the j value;
judging whether j is less than or equal to the length Nl of the list;
if j is less than or equal to the length Nl of the list, returning to the step of utilizing the BlockMerge function according to the formula U j ,∑ j ,V j ＝BlockMerge(l Ut (j),l ∑t (j),l Vt (j),l Ut (j+1),l ∑t (j+1),l Vt (j + 1)), U is obtained j ,∑ j ,V j ”；
If j is larger than the length Nl of the list, judging whether the length Nl of the list is an odd number;
if the list length Nl is odd, the list l is added Ut Is added to the list l U Will listl ∑t Is added to the list l ∑ Will list l Vt Is added to the list l V ；
If the list length Nl is an even number, adding 1 to the value of i;
judging whether i is less than or equal to the iteration layer level;
if i is less than or equal to the level of the iteration layer, returning to the established empty list l Ut 、l ∑t And l Vt Use list l U Is 1 of Ut Assign value, list l ∑ Is 1 of ∑t Assign value, list l V Is 1 Vt Assigning values ";
if i is greater than the iteration layer level, ending the loop;
will list l U 、l ∑ And l V Output as a first singular value decomposition result;
performing singular value decomposition on the aggregation result to obtain a second singular value decomposition result, specifically including:
according to formula N c = round (n/col' + 0.45) get cut column block number; n represents the column number of the aggregation result matrix M to be decomposed, and col' represents the width of a column block after preset splitting;
building a list l for storing the cut column blocks A And a list l for storing the decomposition results in each column of blocks U 、l ∑ And l V ；
Cutting matrix M into N columns c Column block, N c Put column blocks into list l A Performing the following steps;
to N c Singular value decomposition is carried out on each column block, and the result is correspondingly stored into a list l U 、l v And l V Performing the following steps;
will l U 、l ∑ And l V The singular value decomposition result of the matrix M is obtained as an input to the domargeofblocks function.
CN202110561327.6A 2021-05-22 2021-05-22 Subspace system identification method and system Active CN113110068B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110561327.6A CN113110068B (en) 2021-05-22 2021-05-22 Subspace system identification method and system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110561327.6A CN113110068B (en) 2021-05-22 2021-05-22 Subspace system identification method and system

Publications (2)

Publication Number Publication Date
CN113110068A CN113110068A (en) 2021-07-13
CN113110068B true CN113110068B (en) 2023-03-10

Family

ID=76722965

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110561327.6A Active CN113110068B (en) 2021-05-22 2021-05-22 Subspace system identification method and system

Country Status (1)

Country Link
CN (1) CN113110068B (en)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107632522A (en) * 2017-08-31 2018-01-26 南京理工大学 One proton exchanging film fuel battery Nonlinear state space model discrimination method
CN112684707A (en) * 2020-12-25 2021-04-20 扬州大学 Styrene bulk polymerization anti-interference distribution shape control method based on interference observer

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8346712B2 (en) * 2009-11-24 2013-01-01 King Fahd University Of Petroleum And Minerals Method for identifying hammerstein models
US9381823B2 (en) * 2014-07-17 2016-07-05 Ford Global Technologies, Llc Real-time battery estimation

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107632522A (en) * 2017-08-31 2018-01-26 南京理工大学 One proton exchanging film fuel battery Nonlinear state space model discrimination method
CN112684707A (en) * 2020-12-25 2021-04-20 扬州大学 Styrene bulk polymerization anti-interference distribution shape control method based on interference observer

Also Published As

Publication number Publication date
CN113110068A (en) 2021-07-13

Similar Documents

Publication Publication Date Title
KR20190073303A (en) Method and electronic device for convolution calculation in neutral network
KR20190073302A (en) Method and electronic device for convolution calculation in neutral network
JP2019537139A5 (en)
Harrison et al. Passage time distributions in large Markov chains
CN109033021B (en) Design method of linear equation solver based on variable parameter convergence neural network
Ahmed et al. An assessment of some solvers for saddle point problems emerging from the incompressible Navier–Stokes equations
Ren et al. AIPerf: Automated machine learning as an AI-HPC benchmark
KR20190093932A (en) Arithmetic processing apparatus and method in deep running system
CN111640296A (en) Traffic flow prediction method, system, storage medium and terminal
Saura et al. Computational kinematics of multibody systems: Two formulations for a modular approach based on natural coordinates
CN113110068B (en) Subspace system identification method and system
CN109540089A (en) It is a kind of based on Bayes-Kriging model bridge elevation approximating method
Khimich et al. Numerical study of the stability of composite materials on computers of hybrid architecture
Skalicky et al. Performance modeling of pipelined linear algebra architectures on FPGAs
Bahalkeh et al. Efficient system matrix calculation for manufacturing systems
CN114462524A (en) Clustering method for data center batch processing operation
Mohammadi Parallel reverse time integration and reduced order models
Vabishchevich NUMERICAL SOLUTION OF NONSTATIONARY PROBLEMS FOR A CONVECTION AND A SPACE-FRACTIONAL DIFFUSION EQUATION.
Hajgató et al. Accelerating convergence of fluid dynamics simulations with convolutional neural networks
Sima et al. Numerical Investigation of Newton's Method for Solving Discrete-time Algebraic Riccati Equations.
US20200371838A1 (en) Scheduling operations
CN114691345A (en) Calculation framework suitable for SLAM nonlinear parallelization chip and working method
Varga Solution of positive periodic discrete Lyapunov equations with applications to the balancing of periodic systems
Kuźnik et al. Grammar-Based Multi-Frontal Solver for One Dimensional Isogeometric Analysis with Multiple Right-Hand-Sides
Skalicky et al. Performance modeling of pipelined linear algebra architectures on FPGAs

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant