DATA PROCESSING EQUIPMENT, CONTROL SYSTEM, DATA PROCESSING METHOD AND PROGRAM

Information

  • Patent Application
  • 20240403384
  • Publication Number
    20240403384
  • Date Filed
    September 02, 2022
    2 years ago
  • Date Published
    December 05, 2024
    17 days ago
Abstract
A data processing device to estimate the limit:
Description
TECHNICAL FIELD

This disclosure relates to a data processing device, a control system, a data processing method and a program.


BACKGROUND ART

A controllability Gramian is known as an indicator for controllability of a system, see, for example, Non-Patent Document 1-6. Whether a system is controllable can be determined by whether the controllability Gramian is regular. Furthermore, the magnitude of the eigenvalues of the controllability Gramian quantitatively indicates the degree of influence of the input on the state.


PRIOR ART LITERATURE
Non-Patent Literature



  • Non-Patent Literature 1: F. Pasqualetti, S. Zampieri, and F. Bullo, “Controllability metrics, limitations and algorithms for complex networks,” IEEE Transactions on Control of Network Systems, vol. 1, no. 1, pp. 40-52, 2014.

  • Non-Patent 2: V. Tzoumas, M. A. Rahimian, G. J. Pappas, and A. Jadbabaie, “Minimal actuator placement with bounds on control effort,” IEEE Transactions on Control of Network Systems, vol. 3, no. 1, pp. 40-52, 2016.

  • Non-patent 3: K. Kashima, “Noise response data reveal novel controllability Gramian for nonlinear network dynamics,” Scientific Reports vol. 6, Art. 27300, 2016.

  • Non-Patent 4: V. M. Preciado and M. A. Rahimian, “Controllability Gramian spectra of random networks,” 2016 American Control Conference, pp. 3874-3879, 2016.

  • Non-Patent 5: X. Cheng and J. M. A. Scherpen, “A new controllability Gramian for semistable systems and its application to approximation of directed networks,” 56th Annual Conference on Decision and Control, pp. 3823-3828, 2017.

  • Non-Patent 6: S. Zhao and F. Pasqualetti, “Networks with diagonal controllability Gramians: analysis, graphical conditions, and design algorithms,” Automatica, vol. 102, pp. 10-18, 2019.

  • Non-Patent 7: Z. Wang and D. Liu, “Data-based controllability and observability analysis of linear discrete-time systems,” IEEE Transactions on Neural Networks, vol. 22, no. 12, pp. 2388-2392, 2011.

  • Non-Patent 8: D. Bhattacharjee, B. Klose, G. B. Jacobs, and M. S. Hemati, “Data-driven selection of actuators for optimal control of airfoil separation,” Theoretical and Computational Fluid Dynamics, vol. 34, no. 4, pp. 557-575, 2020.

  • Non-Patent Document 9: K. Zhou, J. C. Doyle, and K. Glover, Robust and Optimal Control, Prentice Hall, 1996.

  • Non-Patent Document 10: C. Kenney and G. Hewer, “Trace norm bounds for stable Lyapunov operators,” Linear Algebra and its Applications, vol. 221, pp. 1-18, 1995.



SUMMARY OF INVENTION
Technical Problem

A Controllability Gramian can be easily calculated if a mathematical model of a system is known. Hereinafter, calculation methods based on the mathematical model of a system are called “model-based methods”. However, the amount of data necessary for modeling the system is not always available. In such cases, a mathematical model of the system cannot be obtained.


In contrast to the model-based method, calculation methods that calculate based on data of a state trajectory of a system without using a mathematical model are called “data-driven methods”. Data-driven methods have the advantage of requiring fewer decision variables than model-based methods because they do not need to identify a mathematical model of the system.


Methods for estimating controllability Gramians, see, for example, Non-Patent Document 7, and maximizing controllability Gramians, see, for example, Non-Patent Document 8, using data-driven methods are known. Both of these methods are discrete-time models formulated on the basis of discrete time. However, the physical information contained in the data is represented clearly in the continuous-time model, which is formulated on the basis of continuous time, whereas the discrete-time model has a problem of unclarity. In other words, conventional discrete-time models have the problem that it is difficult to utilize prior knowledge about the characteristics of continuous-time systems.


The purpose of this disclosure is to estimate the controllability Gramian for a system of which mathematical model is unknown, using a data-driven method for continuous-time systems.


Solution to Problem

In order to solve the above problem, a data processing device according to an embodiment of this disclosure estimates the limit:









[

equation


140

]









G

(

)










of a controllability Gramian:









[

equation


138

]









G

(
t
)










defined by:









[

equation


2

]











G

(
t
)

:

=



0


t




?


BB



?


















?

indicates text missing or illegible when filed




in:









[

equation


144

]









t
=











when:









[

equation


1

]











x
.

(
t
)

=


Ax

(
t
)

+

Bu

(
t
)












holds, where:









[

equation


142

]









x

(
t
)










is an n-dimensional vector representing the state of a control object:









[

equation


139

]









u

(
t
)










is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix.


This data processing device comprises:

    • a data acquisition unit that acquires a set of time-series state data:









[

equation


5

]










x

(


[


t

1

1


,

t

1

2



]

,

x

1

1



)

,

x

(


[


t

2

1


,

t

2

2



]

,

x

1

1



)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)











for the following q time intervals:









[

equation


4

]










[


t

i

1


,

t

i

2



]




(


i
=
1

,
2
,


,
q

)











when:









[

equation


39

]










u

(
t
)


0











holds;

    • a controllability Gramian calculation unit that defines:









[

equation


75

]










z

(
t
)



R
n











expressed as:









[

equation


74

]











z
˙

(
t
)

=


A




z

(
t
)






(
13
)







based on the sets of state data:









[

equation


145

]









E

(
t
)










acquired by the data acquisition unit 10, calculates:









[

equation


81

]













z


(

t
2

)



Xz

(

t
2

)


-



z


(

t
1

)



Xz

(

t
1

)



=

-




t
1




t
2






z


(
t
)



BB




z

(
t
)


dt







(
16
)












[

equation


89

]










z

(


t
+

t

i

1



,

t

i

1


,

x

i

1



)

=



(


E

(
t
)



E
0

-
1



)





x

i

1













and estimates:









[

equation


146

]










G

(

)

=
X










by numerically obtaining the solution X of the following linear equation:









[

equation


6

]






















x

i

1



(


E

(
h
)



E
0

-
1



)




X

(


E

(
h
)



E
0

-
1



)





x

i

1



-


x

i

1





Xx

i

1




=

-



0


h





x

i

1



(


E

(

?

)



E
0

-
1



)





BB


(


E

(

?

)



E
0

-
1



)





x

i

1



dt









(


i
=
1

,
2
,


,
q

)







?

indicates text missing or illegible when filed




with respect to:









[

equation


157

]



















E

(
t
)

:=



[




x


(



?

+

?


,

?

,

x
11


)





x


(



?

+


?

21


,


?

21

,

x
21


)









x


(


t
+

t

n

1



,

t

n

1


,

x

n

1



)





]











[

equation


158

]











E
0

:=

[




x
11




x
21







x

n

1





]


;













?

indicates text missing or illegible when filed




and

    • an output unit that outputs the controllability Gramian estimated by the controllability Gramian calculation unit.


In one embodiment, the set of time-series state data acquired by the data acquisition unit may include noise and the controllability Gramian calculation unit estimates:









[

equation


146

]










G

(

)

=
X










by numerically obtaining the solution X of the following liner equation:









[

equation


108

]











(
23
)













?


(


E

(
h
)



E
0
+


)




X

(


E

(
h
)



E
0
+


)





x

i

1



-


x

i

1




X

?



=

-



0


h





x

i

1



(


E

(
t
)



E
0
+


)





BB


(


E

(
t
)



E
0
+


)





x

i

1



dt









(


i
=
1

,
2
,


,
p

)







?

indicates text missing or illegible when filed




instead of the liner equation:









[

equation


6

]






















x

i

1



(


E

(
h
)



E
0

-
1



)




X

(


E

(
h
)



E
0

-
1



)





x

i

1



-


x

i

1





Xx

i

1




=

-



0


h





x

i

1



(


E

(
t
)



E
0

-
1



)





BB


(


E

(
t
)



E
0

-
1



)





x

i

1



dt










(


i
=
1

,
2
,


,
q

)

.




In one embodiment, the controllability Gramian calculation unit performs numerical calculations using prior knowledge about the signs of some or all of the matrix components of the solution X.


Another embodiment of this disclosure is also a data processing device. This data processing device estimates the matrix B which maximizes the trace:









[

equation


141

]









tr

(

G

(

)

)










of the limit:









[

equation


140

]









G

(

)










of a controllability Gramian:









[

equation


138

]









G

(
t
)










defined by:









[

equation


2

]










G

(
t
)

:=



0


t




?


BB



?

d

τ















?

indicates text missing or illegible when filed




in:









[

equation


144

]









t
=











when:









[

equation


1

]











x
˙

(
t
)

=


Ax

(
t
)

+

Bu

(
t
)












holds, where:









[

equation


142

]









x

(
t
)










is an n-dimensional vector representing the state of the control object:









u

(
t
)




[

equation


139

]







is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix. This data processing device comprises:

    • a data acquisition unit that acquires a set of time-series state data:










x

(


[


t

1

1


,

t

1

2



]

,

x

1

1



)

,

x

(


[


t

2

1


,

t

2

2



]

,

x

1

1



)

,


,

x

(


[


t

q

1


,

t

q

2



]

,


x

q

1



)





[

equation


5

]







for the following q time intervals:










[


t

i

1


,

t

i

2



]




(


i
=

1

,
2
,


,
q

)





[

equation


4

]







when:










u

(
t
)


0




[

equation


39

]







holds;

    • a maximization condition calculation unit that estimates the matrix B which maximizes:









t


r

(

G

(

)

)





[

equation


141

]







by numerically obtaining the solution:










Y
~



*





[

equation


9

]







of the following linear equation:














x


(


t

i

2


,

t

i

1


,

x

i

1



)



Y


x

(


t

i

2


,

t

i

1


,

x

i

1



)


-


x

i

1




Y


x

i

1




=

-




t

i

1





t

i

2







x


(

t
,

t

i

1


,

x

i

1



)



Q


x

(

t
,

t

i

1


,

x

i

1



)



dt








(


i
=
1

,
2
,


,
q

)





[

equation


8

]







when:









Q
:=
I




[

equation


7

]







holds;


and

    • an output unit that outputs the input matrix when the controllability Gramian is maximized based on the estimated maximization condition.


In one embodiment, the maximization condition calculation unit may obtain a first input matrix:










B
˜

1


*





[

equation


11

]







by calculating the unit eigenvector corresponding to the maximum eigenvalue of:










Y
˜



*





[

equation


9

]







when:











B




m
=
1







=

R

n
×
m







[

equation


47

]







holds.


In one embodiment, the maximization condition calculation unit may obtain a second input matrix:










B
˜

2


*





[

equation


12

]







by calculating the n×n matrix in which the (k, k) component is 1 and the other components are 0, where the (k, k) component is the one that is the maximum among the diagonal components of:











Y
~



*


,




[

equation


9

]







when:









B
=

D
n





[

equation


49

]







holds.


Yet another embodiment of this disclosure is a control system. This control system controls external control objects. This control system comprises:

    • a sensor that detects a set of time-series state data from the control objects;
    • the aforementioned data processing device; and
    • a control unit to control the control objects,
    • wherein the sensor transmits the detected set of state data to the data acquisition unit of the data processing device,
    • wherein the data processing device transmits the estimated:










G

(

)

=
X




[

equation


146

]









    • to the control unit and

    • wherein the control unit controls the control objects based on:













G

(

)

=

X
.





[

equation


146

]







Yet another embodiment of this disclosure is a data processing method. This method estimates the limit:









G

(

)




[

equation


140

]







of a controllability Gramian:









G

(
t
)




[

equation


138

]







defined by:










G

(
t
)

:=



0
t



e

A

τ




BB




e


A



τ



d

τ






[

equation


2

]







in:









t
=





[

equation


144

]







when:











x
˙

(
t
)

=


A


x

(
t
)


+

B


u

(
t
)







[

equation


1

]







holds, where:









x

(
t
)




[

equation


142

]







is an n-dimensional vector representing the state of a control object:









u

(
t
)




[

equation


139

]







is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix. This method comprises:

    • acquiring a set of time-series state data:










x

(


[


t

1

1


,

t

1

2



]

,

x

1

1



)

,

x

(


[


t
21

,

t

2

2



]

,

x

1

1



)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)





[

equation


5

]







for the following q time intervals:










[


t

i

1


,

t

i

2



]




(


i
=
1

,
2
,


,
q

)





[

equation


4

]







when:










u

(
t
)


0




[

equation


39

]







holds;

    • defining:










z

(
t
)



R
n





[

equation


75

]







expressed as:











z
˙

(
t
)

=


A




z

(
t
)






[

equation


74

]







based on the sets of state data:









E

(
t
)




[

equation


145

]







acquired by the data acquisition unit 10, calculating:









[

equation


81

]













z


(

t
2

)



Xz

(

t
2

)


-



z


(

t
1

)



Xz

(

t
1

)



=

-




t
1


t
2





z


(
t
)



BB




z

(
t
)


dt







(
16
)













z

(


t
+

t

i

1



,

t

i

1


,

x

i

1



)

=



(


E

(
t
)



E
0

-
1



)





x

i

1







[

equation


89

]







and estimating:










G

(

)

=
X




[

equation


146

]







by numerically obtaining the solution X of the following linear equation:













x

i

1



(


E

(
h
)



E
0

-
1



)




X

(


E

(
h
)



E
0

-
1



)





x

i

1



-


x

i

1




X



x

i

1




=

-





(
)

h




x

i

1



(


E

(
t
)



E
0

-
1



)





BB


(


E

(
t
)



E
0

-
1



)





x

i

1




dt



(


i
=
1

,
2
,

,
q

)








[

equation


6

]







with respect to:










E

(
t
)

:=

[


x

(


t
+

t

1

1



,

t
11

,

x

1

1



)




x

(


t
+

t
21


,

t
21

,

x
21


)







x

(


t
+

t

n

1



,

t



n

1



,

x

n

1



)


]





[

equation


157

]














E
0

:=

[


x
11




x
21







x

n

1



]


;




[

equation


158

]







and

    • outputting the estimated controllability Gramian.


Yet another embodiment of this disclosure is also a data processing method. This method estimates the matrix B which maximizes the trace:









tr

(

G

(

)

)




[

equation


141

]







of the limit:









G

(

)




[

equation


140

]







of a controllability Gramian:









G

(
G
)




[

equation


138

]







defined by:










G

(
t
)

:=



0
t



e

A
τ




BB




e


A



τ



d

τ






[

equation


2

]







in:









t
=





[

equation


144

]







when:











x
˙

(
t
)

=


Ax

(
t
)

+

Bu

(
t
)






[

equation


1

]







holds, where:









x

(
t
)




[

equation


142

]







is an n-dimensional vector representing the state of the control object:









u

(
t
)




[

equation


139

]







is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix. This method comprises:

    • acquiring a set of time-series state data:










x

(


[


t

1

1


,

t

1

2



]

,

x

1

1



)

,

x

(


[


t

2

1


,

t

2

2



]

,

x

1

1



)

,

,

x

(


[


t
q1

,

t

q

2



]

,

x

q

1



)





[

equation


5

]







for the following q time intervals:










[


t

i

1


,

t

i

2



]




(


i
=
1

,
2
,

,
q

)





[

equation


4

]







when:










u

(
t
)


0




[

equation


39

]







holds;

    • estimating the matrix B which maximizes:









tr

(

G

(

)

)




[

equation


141

]







by numerically obtaining the solution:






[

equation


9

]







Y
˜

*




of the following linear equation:






[

equation


8

]










x


(


t

i

2


,

t

i

1


,

x

i

1



)



Yx

(


t

i

2


,

t

i

1


,

x

i

1



)


-


x

i

1





Yx

i

1




=

-




t

i

1





t

i

2







x



(

t
,

t

i

1


,

x

i

1



)



Qx

(

t
,

t

i

1


,

x

i

1



)


dt









(


i
=
1

,
2
,


,
q

)




when:






[

equation


7

]






Q
:=
I




holds;


and

    • outputting the input matrix when the controllability Gramian is maximized based on the estimated maximization condition.


Another embodiment of this disclosure is a program. This program estimates the limit:






[

equation


140

]






G

(

)




of a controllability Gramian:






[

equation


138

]






G

(
t
)




defined by:






[

equation


2

]







G

(
t
)

:=



0


t




e

A

τ




BB




e

A


τ




d

τ






in:






[

equation


144

]






t
=





when:






[

equation


1

]








x
˙

(
t
)

=



A
x

(
t
)

+

Bu

(
t
)






holds, where:






[

equation



1

42


]






x

(
t
)




is an n-dimensional vector representing the state of a control object:






[

equation


139

]






u

(
t
)




is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix. This program causes the computer to perform the method, comprising:


acquiring a set of time-series state data:






[

equation


5

]







x

(


[


t
11

,

t

1

2



]

,

x

1

1



)

,

x

(


[


t

2

1


,

t

2

2



]

,

x

1

1



)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)





for the following q time intervals:






[

equation


4

]







[


t

i

1


,

t

i

2



]




(


i

=
1

,
2
,


,
q

)





when:






[

equation


39

]







u

(
t
)


0




holds;


defining:






[

equation


75

]







z

(
t
)




R
n





expressed as:






[

equation


74

]











z
˙

(
t
)

=


A




z

(
t
)







(
13
)








based on the sets of state data:






[

equation


145

]






E

(
t
)




acquired by the data acquisition unit 10, calculating:






[

equation






81

]













z


(

t
2

)


X


z

(

t
2

)


-



z


(

t
1

)


X


z

(

t
1

)



=

-




t
1




t
2






z


(
t
)



BB




z

(
t
)


dt








(
16
)










[

equation


89

]







z

(


t
+

t

i

1



,

t

i

1


,

x

i

1



)

=



(


E

(
t
)



E
0

-
1



)





x

i

1







and estimating:






[

equation


146

]







G

(

)

=
X




by numerically obtaining the solution X of the following linear equation:






[

equation


6

]










x

i

1



(


E

(
h
)



E
0

-
1



)




X

(


E

(
h
)



E
0

-
1



)





x

i

1




-


x

í


1






Xx

i

1




=


-




0


h





x

í

1



(


E

(
t
)



E
0

-
1



)





BB


(


E

(
t
)



E
0

-
1



)





x

i

1



dt









(


i
=
1

,
2
,


,
q

)




with respect to:






[

equation


157

]







E

(
t
)

:=

[


x

(


t
+

t
11


,

t
11

,

x
11


)




x

(


t
+

t
21


,

t
21

,

x
21


)







x

(


t
+

t

n

1



,

t

n

1


,

x

n

1



)


]







[

equation


158

]








E
0

:=

[


x
11



x
21







x

n

1



]


;




and


outputting the estimated controllability Gramian.


Yet another embodiment of this disclosure is also a program. This program estimates the matrix B which maximizes the trace:






[

equation


141

]






t


r

(

G

(

)

)





of the limit:






[

equation


140

]






G

(

)




of a controllability Gramian:






[

equation


138

]






G

(
t
)




defined by:






[

equation


2

]







G

(
t
)

:=



0


t




e

A

τ




BB




e


A



τ



d

τ






in:






[

equation


144

]






t
=





when:






[

equation


1

]








x
˙

(
t
)

=


Ax

(
t
)

+

Bu

(
t
)






holds, where:






[

equation


142

]






x

(
t
)




is an n-dimensional vector representing the state of the control object:






[

equation


139

]






u

(
t
)




is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix. This program causes the computer to perform the method, comprising:

    • acquiring a set of time-series state data:






[

equation


5

]







x

(


[


t

1

1


,

t

1

2



]

,

x

1

1



)

,

x

(


[


t

2

1


,

t

2

2



]

,

x

1

1



)

,


,

x

(


[


t

q

1


,


t

q

2



]

,

x

q

1



)





for the following q time intervals:






[

equation


4

]







[


t

i

1


,

t

i

2



]




(


i
=
1

,
2
,


,
q

)





when:






[

equation


39

]







u

(
t
)


0




holds;

    • estimating the matrix B which maximizes:






[

equation


141

]






tr

(

G

(

)

)




by numerically obtaining the solution:






[

equation


9

]







Y
˜

*




of the following linear equation:






[

equation


8

]










x


(


t

i

2


,

t

i

1


,

x

i

1



)



Yx

(


t

i

2


,

t

i

1


,

x

i

1



)


-


x

i

1




Y


x

i

1





=


-
⁠⁠




t

i

1





t

i

2







x



(

t
,

t

i

1


,

x

i

1



)



Qx


(

t
,

t

i

1


,

x

i

1



)



dt









(


i
=
1

,
2
,


,
q

)




when:






[

equation


7

]






Q
:=
I




holds;


and

    • outputting the input matrix when the controllability Gramian is maximized based on the estimated maximization condition.


Yet another embodiment of this disclosure is a data processing device. This data processing device estimates:






[

equation


232

]







G
c

(

A
+
Δ

)




when:






[

equation


1

]








x
.

(
t
)

=


Ax

(
t
)

+

Bu

(
t
)






holds and the limit of a controllability Gramian of a matrix A, which is defined by:






[

equation


2

]








G

(
t
)

:=



0


t




e

A

τ




BB




c


A



τ



d

τ



,




is defined by:






[

equation


231

]








G
c

(
A
)

,




where:






[

equation


142

]






x

(
t
)




is an n-dimensional vector representing the state of a control object:






[

equation


139

]






u

(
t
)




is an m-dimensional vector representing the control input, A is an unknown n×n matrix, is a known n×m matrix and 4 the amount of change in A. This data processing device comprises:

    • a data acquisition unit that acquires a set of time-series state data:






[

equation


5

]







x

(


[


t
11

,

t

1

2



]

,

x
11


)

,

x

(


[


t
21

,

t

2

2



]

,

x

1

1



)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)





for the following q time intervals:






[

equation


4

]







[


t

i

1


,

t

i

2



]




(


i
=
1

,
2
,


,
q

)





when:






[

equation


39

]







u

(
t
)


0




holds;

    • a controllability Gramian calculation unit that defines:






[

equation


75

]







z

(
t
)




R
n





expressed as:











z
˙

(
t
)

=


A




z

(
t
)






[

equation


74

]







based on the sets of state data:









E

(
t
)




[

equation


145

]







acquired by the data acquisition unit 10, calculates:









[

equation


81

]













z


(

t
2

)



Xz

(

t
2

)


-



z


(

t
1

)



Xz

(

t
1

)



=

-




t
1


t
2





z


(
t
)



BB




z

(
t
)


dt







(
16
)













z

(


t
+

t

i

1



,

t

i

1


,

x

i

1



)

=



(


E

(
t
)



E
0

-
1



)





x

i

1







[

equation


89

]







and estimates:











G
c

(

A
+
Δ

)

=
X




[

equation


233

]







by numerically obtaining the solution X of the following linear equation.









[

equation


225

]













z
i


(
h
)


X



z
i

(
h
)


-



z
i


(
0
)




Xz
i

(
0
)


+



0
h




z
i


(
t
)



(


Δ

X

+

X


Δ




)




z
i

(
t
)


dt



=

-



0
h




z
i


(
t
)



BBz

(
t
)


dt



(


i
=
1

,
2
,


,
N

)








(
47
)







with respect to:










E

(
t
)

:=

[







x


(

t
+










t
11

,

t
11

,

x
11


)










x


(

t
+










t
11

,

t
11

,

x
11


)













x


(

t
+










t
11

,

t
11

,

x
11


)







]





[

equation


157

]














E
0

:=

[










x
11




x
21













x

n

1





]


;




[

equation


158

]







and


an output unit that outputs the estimated controllability Gramian.


Any combination of the above components, and any transformation of the expression of the invention among devices, methods, systems, recording media, computer programs, etc., is also valid as a form of the invention.


Advantageous Effects of Invention

According to the present disclosure, the controllability Gramian for a system of which mathematical model is unknown can be estimated, using a data-driven method for continuous-time systems.





BRIEF DESCRIPTION OF DRAWINGS


FIG. 1 is a functional block diagram of a data processing system according to the first embodiment;



FIG. 2 is a functional block diagram of a data processing system according to the second embodiment;



FIG. 3 is a functional block diagram of a control system according to the seventh embodiment;



FIG. 4 is a functional block diagram of a control system according to the eighth embodiment;



FIG. 5 is a flowchart showing the processing procedure of a data processing method according to the ninth embodiment;



FIG. 6 is a flowchart showing the processing procedure of a data processing method according to the tenth embodiment;



FIG. 7 shows time series data without noise in the time interval:










[


t
11

,

t

1

2



]

;




[

equation


151

]








FIG. 8 shows time series data without noise in the time interval:










[


t
21

,

t

2

2



]

;




[

equation


152

]








FIG. 9 shows time series data without noise in the time interval:










[


t
31

,

t

3

2



]

;




[

equation


153

]








FIG. 10 shows time series data without noise in the time interval:










[


t

4

1


,

t

4

2



]

;




[

equation


154

]








FIG. 11 shows time series data without noise in the time interval:










[


t

5

1


,

t

5

2



]

;




[

equation


155

]








FIG. 12 shows time series data without noise in the time interval:










[


t
61

,

t

6

2



]

;




[

equation


156

]








FIG. 13 shows time series data with noise in the time interval:










[


t
11

,

t
12


]

;




[

equation


151

]








FIG. 14 shows time series data with noise in the time interval:










[


t

2

1


,

t

2

2



]

;




[

equation


152

]








FIG. 15 shows time series data with noise in the time interval:










[


t

3

1


,

t

3

2



]

;




[

equation


153

]








FIG. 16 shows time series data with noise in the time interval:










[


t
41

,

t

4

2



]

;




[

equation


154

]








FIG. 17 shows time series data with noise in the time interval:










[


t
51

,

t

5

2



]

;




[

equation


155

]








FIG. 18 shows time series data with noise in the time interval:










[


t
61

,

t

6

2



]

;




[

equation


156

]







and



FIG. 19 is a functional block diagram of a control system for a variant according to the seventh embodiment.





DESCRIPTION OF EMBODIMENTS

The disclosure will be explained below with reference to drawings based on suitable embodiments. The embodiments are examples rather than limitations of the disclosure. All features or combinations of features described in the embodiments are not necessarily essential to the disclosure. Identical or equivalent components, parts, and processes shown in each drawing shall be given the same symbol, and redundant explanations will be omitted where appropriate. The scale and shape of each part shown in each drawing are set for convenience in order to facilitate explanation, and are not to be construed as limiting unless otherwise noted. When terms such as “first,” “second,” etc. are used in this specification or in the claims, unless otherwise mentioned, these terms do not indicate any order or degree of importance, but are intended only to distinguish one configuration from another. In addition, in each drawing, some parts of the components that are not important in explaining the embodiments are omitted.


Before describing the specific embodiments, let us first explain the basic findings.


Definition of Symbols

The following are definitions of mathematical symbols used herein.


R: real number.










R
+

:




[

equation


128

]







positive real number.










R

0
+


:




[

equation


129

]







non-negative real number.










S
n

:




[

equation


130

]







n×n symmetric matrix.










D
n

:




[

equation


131

]







n×n diagonal matrix.


0: zero matrix.


I: unit matrix.










e
i

:




[

equation


132

]







the i-th standard basis of the n-dimensional real vector space.










A
+

:




[

equation


133

]







Moore-Penrose pseudo-inverse of matrix A.











λ
max

(
A
)

:




[

equation


134

]







the maximum eigenvalue of real symmetric matrix of which all eigenvalues are real numbers.










tr

(
A
)

:




[

equation


135

]







trace of square matrix A.


Note that for any square matrices A and B:










tr

(
AB
)

=

tr

(
BA
)





[

equation


13

]







holds.












A


F

:




[

equation


14

]







Frobenius norm of matrix A, i.e.:












A


F

=


tr

(


A



A

)






[

equation


15

]







holds.


Furthermore, for an n×n symmetric matrix P of which (i, j) components are:










p
ij

,




[

equation


136

]







i.e.:











p
ij

=

p
ji


,




[

equation


137

]







column vector:










rvec

(
P
)



R


(

1
/
2

)



n

(

n
+
1

)







[

equation


16

]







and row vector:










lvec

(
x
)



R


(

1
/
2

)



n

(

n
+
1

)







[

equation


17

]







are defined as follows:










rvec

(
P
)

=

[


p
11





"\[LeftBracketingBar]"





p
12




p
22






"\[RightBracketingBar]"




|




p

1

n





p

2

n










p
nn

]












[

equation


18

]














x



Px

=


lvec

(
x
)




rvec

(
P
)

.






[

equation


19

]







For example:










lvec

(
x
)

=

[




x
1
2




2


x
1



x
2





x
2
2




]





[

equation


22

]













rvec

(
P
)

=


[




p

1

1





p
12




p
22




]







[

equation


23

]







holds for:









x
=


[




x
1






x
2




]



R
2






[

equation


20

]












P
=


[




p

1

1





p
12






p
21




p
22




]



S
2






[

equation


21

]







respectively.


[Linear Systems and controllability Gramians]


We consider the following linear system.






[

equation


24

]












x
˙

(
t
)

=


Ax

(
t
)

+

Bu

(
t
)








(
1
)








Here:









x

(
t
)



R
n





[

equation


25

]







is the state of the system:










u

(
t
)



R
m





[

equation


26

]







is the control input:









A


R

n
×
n






[

equation


27

]







is an n×n constant matrix and:









B


R

n
×
m






[

equation


28

]







is an n×m constant matrix.


In this case, the controllability Gramian:









G

(
t
)




[

equation


138

]







is defined as follows.










G

(
t
)

:=



0
t



e

A
τ




BB




e


A



τ



d

τ






[

equation


29

]







Here:








t


R
+





[

equation


30

]







holds.


The controllability Gramian represents the set of states that can be reached at time t, starting from the initial state given a control input:










u
(
t

)




[

equation


139

]







satisfying the following.






[

equation


31

]








0


t







u

(
τ
)



2

d

τ

=
1









[

equation


138

]







G

(
t
)




is regular when the system is controllable, i.e.:






[

equation


32

]







rank
(

[

B


AB






A

n
-
1



B

]

)

=
n




holds.


The set of states that can be reached in this case is expressed as below.






[

equation


33

]






{

x





R
n





"\[LeftBracketingBar]"




x





G

-
1


(
t
)


x



<
¯


1




}




In particular, the direction that can be reached can be quantified by the limit of controllability Gramians:






[

equation


140

]






G

(

)




in the following.






[

equation






140

]






t
=








[

equation


141

]







tr

(

G

(

)

)




is the sum of the eigenvalues of matrix below.






[

equation






140

]






G

(

)




Therefore:





[

equation






141

]






t


r

(

G

(

)

)





is a measure of the strength of controllability of the system.


The following lemma is the basis for the computation of the controllability Gramian, see, for example, Non-Patent Literature 9.


(Lemma 1)

For the system shown in (1), we consider the Lyapunov equation as below.









[

equation


34

]











XA


+
AX

=

-

BB








(
3
)








Here:





[

equation


35

]






X




R

n
×
n






are the unknown variables in equation (1). If matrix A is a Hurwitz matrix, i.e., the real part of all eigenvalues of matrix A is negative, there exists a unique solution X, which is equal to the following.






[

equation


140

]






G

(

)




(End of Lemma 1)
[Data-Driven Estimation and Maximization Problem]

Once the mathematical model of (1) is obtained, the controllability Gramian:






[

equation


140

]






G

(

)




is computable and can be optimized by (2) or (3). However, another approach is needed if a mathematical model is not available. One possible approach is to use measured data on the behavior of the system. This is formulated as follows.


We denote the state (1):






[

equation






142

]






x

(
t
)




by:






[

equation






40

]








x

(

t
,


t
1

,

x
1


)




R
n


,




with respect to:






[

equation


36

]







t
i






R

0
+


(


i

=

1

,
2

)







and
:






[

equation


37

]








x
1





R
n


,




under the conditions of:






[

equation






38

]







x

(

t
1

)

=

x
1







and
:






[

equation






39

]







u

(
t
)


0




Also, we denote the segment of the state trajectory in the time interval under the same boundary conditions and inputs as above:









[


t
1

,

t
2


]




[

equation


41

]







by the following.










x

(


[


t
1

,

t
2


]

,

x
1


)




[


t
1

,

t
2


]

×

R
n






[

equation


42

]







Therefore, the problem of estimating and maximizing:









G

(

)




[

equation


140

]







based on given data can be organized as follows.


(Problem 1)

For the system shown in (1), let matrix A be a Hurwitz matrix and unknown. Let a set of data:










x

(


[


t

i

1


,

t

i

2



]

,

x

i

1



)




(


i
=
1

,
2
,


,
q

)





[

equation


43

]







consisting of segments of q state trajectories is given.


Here:









t

i

1


<

t

i

2






[

equation


44

]







holds.


(i) We estimate:









G

(

)




[

equation


140

]







when a known input matrix B is given.


(ii) We obtain the input matrix:









B

B




[

equation


46

]







that maximizes:









t


r

(

G

(

)

)





[

equation


141

]







under the conditions of the following.












B


F


1




[

equation


45

]







Here:








B





m
=
1




R

n
×
m







[

equation


47

]







is a given pair of possible input matrices.


Note the following four points. First, in Problem 1(i), the input matrix B is assumed to be known. On the other hand, the main objective of Problem 1(ii) is to find the optimal input channel for controlling the system. From this point of view, Problem 1(i) corresponds to a performance analysis of the input matrix B. Namely, in this case, the input matrix B is evaluated to measure the performance of the following.









G

(

)




[

equation


140

]







Second, in general:









G

(

)




[

equation


140

]







increases with the norm of B. Therefore, the condition:












B


F


1




[

equation


45

]







is required for the maximization of:









t


r

(

G

(

)

)





[

equation


141

]







in Problem 1(ii).


Third, although B can be chosen flexibly, the following two settings are essential.









B





m
=
1




R

n
×
m







[

equation


47

]












B
=

D
n





[

equation


49

]







The former is the case where the inputs to the system are completely arbitrary. The latter is useful when the system is a network of n nodes, with one input added to each node. As for the other settings of B, approximate solutions can be obtained from the above settings. This may be obtained, for example, by truncating less influential part of B, as in the low-rank matrix approximation based on singular value decomposition.


Fourth, it is assumed that the resulting set of data is a continuous-time signal and contains no noise. However, even if this assumption does not hold, the method described herein is still applicable to the above problem. This point is discussed below.


[Solving the Data-Driven Lyapunov Equation in Stability Analysis]

We consider the system (1) with B=0.









[

equation


50

]











x
.

(
t
)

=

Ax

(
t
)





(
4
)







In this system, the Lyapunov equation is given by the following.









[

equation


51

]










YA
+


A



Y


=

-
Q





(
5
)







Here:








Q


S
n





[

equation


52

]







is a given symmetric matrix and:









Y


S
n





[

equation


53

]







is the unknown variable in this equation. The unique solution:









Y


S
n





[

equation


53

]







exists if A is a Hurwitz matrix. Note that (5) is different from (3) even though:









Q
=

BB






[

equation


54

]







holds.


We now consider solving (5) when A is unknown and a set of state trajectory data is given. The problem is formulated as follows.


(Problem 2)

We assume that the matrix A is a Hurwitz matrix and is unknown with respect to the system shown in (4). We assume that a symmetric matrix:









Q


S
n





[

equation


52

]







and a set of data consisting of segments of q state orbits:












x


(


[


t

i

1


,

t

i

2



]

,

x

i

1



)





(


i
=
1

,
2
,



;
q


)







[

equation


43

]







is given. Here:










t

i

1


<

t

i

2






[

equation


44

]







holds. We obtain the solution of (5):









Y


R

n
×
n






[

equation


55

]







under this condition.


The solution to Problem 2 is given as follows. We consider the state trajectory x in (4) from the initial state x(0).









[

equation


56

]












x


(
t
)



(

YA
+


A



Y


)



x

(
t
)


=


-


x


(
t
)




Qx

(
t
)






(
6
)







is obtained by making a quadratic form of x(t) with respect to both sides of (5). Furthermore:









[

equation


57

]















t
1





t
2






x


(
t
)



(

YA
+


A



Y


)



x

(
t
)


dt


=

-






t
1





t
2






x


(
t
)



Qx

(
t
)


dt







(
7
)







is obtained by integrating both sides of this equation over the time interval as below.









[


t
1

,

t
2


]




[

equation


41

]







The left side of (7) is equal to:












x


(


t


2

)



Yx

(

t
2

)


-



x


(

t
1

)



Yx

(

t
1

)






[

equation


58

]







because:









[

equation


59

]












x


(
t
)



(

YA
+


A



Y


)



x

(
t
)


=





x


(
t
)


YA

+


A




Yx

(
t
)



=






x


(
t
)


Y



x
.

(
t
)


+




x
.



(
t
)



Yx

(
t
)



=


d
dt



(



x


(
t
)



Px

(
t
)


)








(
8
)







holds from (4). From the above, (7) can be expressed as below.









[

equation


60

]












x


(

t
2

)



Yx

(

t
2

)


=




x


(

t
1

)



Yx

(

t
1

)


=

-




t
1




t
2






x


(
t
)


𝒬


x

(
t
)


dt








(
9
)







The following q equations are obtained by applying the sets of data in Problem 2 to (9).









[

equation


61

]











(
10
)














x


(


t

i

2


,

t

i

1


,

x

i

1



)



Yx

(


t

i

2


,

t

i

1


,

x

i

1



)


-


x

i

1





Yx

i

1




=

-




t

i

1





t

i

2







x


(

t
,

t

i

1


,

x

i

1



)


𝒬


x

(

t
,

t

i

1


,

x

i

1



)


dt









(


i
=
1

,
2
,


,
q

)




Here, (10) can be considered to be a linear equation for an unknown matrix Y because:









[

equation


62

]










x

(


t

i

2


,

t

i

1


,

x

i

1



)

,

x

i

1












and the right side are known constants. In the standard form, (10) is expressed as follows.









[

equation


63

]













D
1



rvec

(
Y
)


=

d
1







(
11
)








Here:








[

equation


64

]












rvec

(
Y
)





R


(

1
/
2

)



n

(

n
+
1

)















are unknown variables. Also:









[

equation


234

]












D
1



R

q
×

(

1
/
2

)



n

(

n
+
1

)














is the matrix given by the following.









[

equation


65

]










D
1

:=

[





lvec

(

x

(


t

12



,

t
11

,

x
11


)

)

-

lvec

(

x
11

)








lvec

(

x

(


t

22



,

t
21

,

x
21


)

)

-

lvec


(

x
21

)














lvcc

(

x


(


t

q

2


,

t

q

1


,

x

q

1



)


)

-

lvcc

(

x

q

1


)





]
















[

equation


66

]












d
1



R


(

1
/
2

)



n

(

n
+
1

)














is a vector formed by arranging the values on the right side of (10) with respect to the following.









[

equation


67

]










i
=
1

,
2
,


,
q










The solution to Problem 2 is thus obtained as follows.


(Lemma 2)

We consider Problem 2. If:









[

equation


68

]









rank



(

D
1

)

=


1
2



n

(

n
+
1

)








(
12
)








holds, there exists a unique solution Y to (10), which is equal to the solution of Problem 2 (End of Lemma 2).


Lemma 2 shows that the solution to Problem 2 can be obtained by solving the linear equation (10), namely (11), using a set of state orbit data for the system shown in (4).


[Estimation of controllability Gramians].


The solution of Problem 1(i) is described below.


[Data Transformation]

According to the aforementioned Lemma 1, Problem 1(i) comes down to solving the Lyapunov equation (3) by a data-driven method. Below, we obtain the solution of (3) based on data using the same method that yielded (10). Similar to (6), we make the quadratic form of the state x(t) when B=0 on the left side of (3) and (1).









[

equation


69

]












x


(
t
)



(


XA


+
AX

)



x

(
t
)


=




x


(
t
)



XA




x

(
t
)


+



x


(
t
)



AXx

(
t
)













However, in this case, unlike in (8), we cannot replace:









[

equation


70

]










A




x

(
t
)













and
:









[

equation


71

]











x


(
t
)


A










by:









[

equation


72

]











x
.

(
t
)

,










in general. This is because x(t) is a state of (1) and:









[

equation


73

]











x
˙

(
t
)




A




x

(
t
)












holds, in general.


Therefore first, we introduce the following system defined by the transpose matrix of matrix A:









[

equation


74

]











z
˙

(
t
)

=


A




z

(
t
)







(
13
)








where:









[

equation


75

]










z

(
t
)



R
n











is the state. In this system, we use the symbol:









[

equation


77

]









z

(

t
,

t
1

,

z
1


)










to express the state:









[

equation


143

]









z

(
t
)










when:









[

equation


76

]










z

(

t
1

)

=


z
1



R
n












holds.


Here, we make the quadratic form:









[

equation


78

]











z


(
t
)



(


XA


+
AX

)



z

(
t
)











which satisfies the following.









[

equation


79

]











(
14
)













z


(
t
)



(


XA


+
AX

)



z

(
t
)


=





z


(
t
)



XA




z

(
t
)


+



z


(
t
)



AXz

(
t
)



=





z


(
t
)


X



z
.

(
t
)


+




z
.



(
t
)



Xz

(
t
)



=


d
dt




z


(
t
)



Xz

(
t
)








On the other hand:









[

equation


80

]














t


1




t
2






z


(
t
)



(


XA


+
AX

)



z

(
t
)


dt


=

-





t


1




t
2






z


(
t
)



BB




z

(
t
)


dt







(
15
)







holds by (3).


Therefore, we obtain:









[

equation


81

]













z


(

t
2

)



Xz

(

t
2

)


-



z


(

t
1

)



Xz

(

t
1

)



=

-





t


1




t
2






z


(
t
)



BB




z

(
t
)


dt







(
16
)







by (14) and (15). This has the same form as (19).


However, in this case, the set of data for the state trajectory in (13) in Problem 1 is not obtained. This difficulty can be resolved by the following result.


(Lemma 3)

We consider:









[

equation


43

]










x

(


[


t

i

1


,

t

i

2



]

,

x

i

1



)




(


i
=

1

,
2
,


,
q

)











in Problem 1. We assume:









[

equation


82

]









q

>
¯

n










holds and choose n arbitrary data from the set of data. Here, n is the dimension of:









[

equation


142

]









x

(
t
)










in (1).


It can be considered that the chosen data are:









[

equation


83

]










x

(


[


t

i

1


,

t

i

2



]

,

x

i

1



)




(


i
=
1

,
2
,


,
n

)











without loss of generality. We assume the following.





[equation 84]






E(t):=[x(t+t11,t11,x11) x(t+t21,t21,x21) . . . x(t+tn1,tn1,xn1)]  (17)









[

equation


85

]










E
0

:=

[




x
11




x
21







x

n

1





]






(
18
)










If
:









[

equation


86

]









E
0










is regular:









[

equation


89

]










z

(


t
+

t

i

1



,

t

i
,
1


,

x

i

1



)

=



(


E

(
t
)



E
0

-
1



)





x

i

1













holds with respect to the following.









[

equation


87

]









i


{

1
,
2
,


,
q

}
















[

equation


88

]









t


[

0
,
h

]











Here:








[

equation


90

]









h
:=


min

i


{

1
,
2
,

,

q

}



(


t

i

2


-

t

i

1



)











holds (End of Lemma 3).


(Proof)

We consider data points:









x

i

1





[

equation


92

]







and state trajectory:









x

(


[


t

j

1


,

t

j

2



]

,


x

j

1



)




[

equation


93

]







with respect to the following which are given.









i




{

1
,
2
,


,
q

}





[

equation


87

]












j




{

1
,
2
,


,
n

}





[

equation


91

]







The product of:










x


(


t
+

t

j

1



,

t

j

1


,


x

j

1



)




[

equation


95

]









and
:









x

i

1





[

equation


92

]







for any:









t




[

0
,



t

j

2


-

t

j

1




]





[

equation


94

]







is given by the following.












x


(


t
+

t

j

1



,

t

j

1


,

x

j

1



)



x

i

1



=




(


e
At



x

j

1



)





x

i

1



=


x

j

1





e

A


t





x

i

1








[

equation


96

]









=


x

j

1





𝓏

(


t
+

t

j

1



,

t

j

1


,

x

i

1



)






By obtaining this relation for j=1, 2, . . . , n, the following is obtained for t.












E


(
t
)



x

i

1



=


E
0




𝓏

(


t
+

t
i1


,

t
i1

,

x

i

1



)






[

equation


97

]












t




[

0
,
h

]





[

equation


88

]







This means that:










𝓏

(


t
+

t

i

1



,

t

i

1


,

x

i

1



)

=



(

E


(
t
)



E
0

-
1



)





x

i

1







[

equation


89

]







holds if:









E
0




[

equation


86

]







is regular (End of proof).


We can transform the set of data in Problem 1 into the set of data in the state trajectory in (13) by this lemma. Note that the assumption of:









q

>
¯

n




[

equation


82

]







is not restrictive. This is because the number of data can be increased by dividing the state trajectory into multiple segments.


[Data-Driven Estimation]

The linear equation:









(
17
)










E

(
t
)

:=

[




x


(


t
+

t
11


,

t
11

,

x
11


)





x


(


t
+

t
21


,

t
21

,

x
21


)








x


(


t
+

t

n

1



,

t

n

1


,

x

n

1



)










[

equation


84

]







is obtained by (16), (19) and the following.












x


(


[


t

i

1


,

t

i

2



]

,

x

i

1



)





(


i
=
1

,
2
,


,
n

)







[

equation


83

]







This is expressed as:









(
21
)











D
2



rvec

(
X
)


=

d
2





[

equation


102

]







with respect to:










d
2





R


(

1
/
2

)



n

(

n
+
1

)







[

equation


101

]







and the following vector.










D
2

:=

[





lvec

(



(

E


(
h
)



E
0

-
1



)





x
11


)

-

lvec

(

x
11

)








lvec

(



(

E


(
h
)



E
0

-
1



)





x
21


)

-

lvec

(

x
21

)













lvec

(



(

E


(
h
)



E
0

-
1



)





x

q

1



)

-

lvec

(

x

q

1


)





]





[

equation


100

]







As a result, the solution to Problem 1 is obtained as follows.


(Theorem 1)

We consider Problem 1(i). We assume the following.











q

n





[

equation


82

]









If
:









[

equation


103

]










rank
(

D
2

)

=


1
2



n

(

n
+
1

)






(
22
)







holds, there exists a unique solution X to (20), which is equal to the solution of Problem 1 (End of Theorem 1).


According to Theorem 1, the solution of Problem 1(i) is given as the solution of the linear equation (20).


In the above discussion, it was assumed that the set of data obtained is a continuous-time signal and contains no noise. However, this is not limited to this, and as shown below, the aforementioned results can also be applied to data containing noise.


We assume that data:












x



(

[


t

i

1


,

t

i

2



]

)


,


x

i

1



)




(


i
=
1

,
2
,



;
q


)





[

equation


43

]







is assumed to contain noise. In order to suppress the effect of noise, it is preferable to use more data to create the matrix:









E

(
t
)




[

equation


145

]









and
:









E
0




[

equation


86

]







in Lemma 3. Therefore:









E

(
t
)



R



n
×
p







[

equation


105

]













E
0



R



n
×
p







[

equation


106

]







are obtained with respect to:










p

<
¯

q

.




[

equation


104

]







Then, from (19) and the properties of the Moore-Penrose pseudo-inverse, the relationship:










z

(


t
+

t

i

1



,

t

i

1


,

x

i

1



)





(


E

(
t
)



E
0
+


)





x

i

1







[

equation


107

]







is obtained. As a result, an approximate solution to Problem 1(i) is obtained from a variant of (20):









[

equation


108

]














x

i

1



(


E

(
h
)



E
0
+


)



X



(


E

(
h
)



E
0
+


)





x

i

1



-


x

i

1




X


x

i

1




=

-



0


h





x

i

1



(


E

(
h
)



E
0
+


)






BB


(


E

(
t
)



E
0
+


)





x

i

1



dt








(


i
=
1

,
2
,


,
p

)






(
23
)








Here:








E

(
t
)




[

equation


145

]









and
:









E
0




[

equation


86

]







are defined with respect to the p data and the following.









h
:=


min

i


{

1
,
2
,



,
p

}



(


t

i

2


-


t

i

1



)





[

equation


109

]







Next, we assume that the sets of data are given as acyclic sample data of state orbitals. The left side of (20) consists of the endpoints of each state trajectory and can be computed from the sample data of the state trajectories. In contrast, the right side of (20) is the integral of the quadratic form of the state orbit. This can be computed approximately from sample data, for example, using the trapezoidal formula. By applying these principles, an approximate solution to the above problem can be obtained.


[Maximizing Controllability Gramians]

The solution of Problem 1(ii) is described below.


[Characterization of Input Matrices]

The following shows the maximization of the controllability Gramian:









G

(

)




[

equation


140

]







for the input matrix B.


(Lemma 4)

For the system shown in (1), we assume that matrix A is a Hurwitz matrix. We consider the following maximization problem:









[

equation


110

]









{





max

B

B




tr

(

G

(

)

)








s
.
t
.




B


F




<
¯


1








(
24
)







where:









B








m
=
1





R



n
×
m








[

equation


47

]







is a given pair of possible input matrices.


We define Y* as follows.









[

equation


111

]











Y
*

:=



0







e


A



τ




e

A

τ



d

τ



,




(
25
)







which is finite positive definite since matrix A is a Hurwitz matrix. The following proposition holds.


(i)


We define the unit eigenvector:










B
1
*



R


n






[

equation


113

]







with respect to the following largest eigenvalue.









Y
*




[

equation


112

]







The matrix:









B
1
*




[

equation


114

]







is the solution of (24) for:









B





m
=
1




R

n
×
m







[

equation


47

]









and










(
26
)










tr

(

G

(

)

)

=



λ
max

(

Y


)



tr

(


B
1





(

B
1


)




)






[

equation


115

]







holds.


(ii)


We assume that:










B
2
*





R

n
×
n






[

equation


116

]







is a diagonal matrix of which (k, k) component is 1 and the other components are 0. Here:









k




{

1
,
2
,


,
n

}





[

equation


117

]







is the row number, also the column number, when the diagonal component of the following matrix.









Y
*




[

equation


112

]












B
2
*




[

equation


118

]







is the solution of (24) for:









B
=

D
n





[

equation


49

]







and is









(
27
)










t


r

(

G

(

)

)


=




max



i



{

1
,
2
,


,
n

}






y

i



(

=

y
k


)






[

equation


119

]







holds. Here:










y
i




R




[

equation


120

]







is the i-th diagonal component of the following.









Y
*




[

equation


112

]







(End of Lemma4)
(Proof)

(i) See, e.g., Non-Patent Document 10.


(ii) We define the i-th diagonal component of matrix B as the following.










b

i




ϵ


R




[

equation


121

]












(
28
)










tr

(

G

(

)

)

=

tr

(



0




e

A

τ





BB
T



e

A


τ




d

τ


)





[

equation


122

]









=




0




tr

(


e

A

τ





BB




e

A


τ




d

τ

)


d

τ


=



0




tr

(


e

A


τ





e

A

τ




BB



)


d

τ









=


tr

(



0




e

A


τ





e

A

τ




BB



d

τ


)

=

tr

(



0




e

A


τ





e

A

τ



d

τ


BB




)






is obtained from (2) and the properties of the trace.


For any:









B




D
n


,
:




[

equation


123

]












BB





[

equation


124

]







is a diagonal matrix. Therefore, from (28):









(
29
)










tr

(

G

(

)

)

=


tr

(



0




e

A


τ





e

A

τ



d

τ


BB




)

=




i
=
1

n




y
i



b
i
2








[

equation


126

]













max



i



{

1
,
2
,


,
n

}





y

i








holds, for any:









B




D
n





[

equation


123

]







which satisfies:












B



F



1




[

equation


45

]









,

i
.
e
.


,
:












[


b
1




b
2







b
n


]




2



<
¯

1.




[

equation


125

]







On the other hand, from the definitions of (29) and:









(
27
)










B
2
*

,


holds



for
:


B
2
*

.







[

equation


118

]







Further:





[

equation


127

]










B
2
*



F

=
1




holds (End of proof).


[Data-Driven Maximization Problem]

Lemma 4 means that, the input matrix that maximizes:






[

equation


141

]






t


r

(

G

(
∞
)

)





is characterized by:






[

equation






112

]






Y
*




in (25). On the other hand:






[

equation


112

]






Y
*




is equal to the unique solution of the Lyapunov equation (5) for:






[

equation


7

]






Q
:=

I
.





Therefore, from Lemma 2 and Lemma 4, the following holds.


(Theorem 2)

We consider Problem 1(ii). We assume that (12) hold with respect to the set of data given in Problem 1. We assume that:






[

equation


9

]







Y
˜

*




is the solution to:






[

equation


61

]













x


(


t

i

2


,

t

i

1


,

x

i

1



)


Y


x

(


t

i

2


,

t

i

1


,

x

i

1



)


-


x

i

1




Y


x

i

1




=




t

i

1



t

i

2






x


(

t
,

t

i

1


,

x

i

1



)


Q


x

(

t
,

t

i

1


,

x

i

1



)


dt



(


i
=
1

,
2
,


,
q

)







(
10
)







with respect to:






[

equation


7

]






Q
:=

I
.





Furthermore, in Lemma 4, we assume that:






[

equation


11

]








B
˜

1
*




and

:








[

equation


112

]






Y
*




are the input matrices with respect to:






[

equation


9

]







Y
~

*




as we defined:






[

equation


118

]






B
2
*




with respect to:






[

equation


114

]







B
1
*




and
:







[

equation


118

]







B
2
*

.




In this case, the following proposition holds.


(i) The input matrix:






[

equation


11

]







B
˜


*




is the solution of Problem 1(ii) with respect to the following:






[

equation


47

]






B





m
=
1




R

n
×
m







(ii) The input matrix:






[

equation


12

]







B
˜

2
*




is the solution of Problem 1(ii) with respect to the following:









B
=

D
n





[

equation


49

]







(End of Theorem 2).

Theorem 2 shows that the solution of Problem 1(ii) can be obtained by solving the linear equation (10) and forming the input matrix by the method of Lemma 4.


The First Embodiment


FIG. 1 is a functional block diagram of a data processing device 1 according to the first embodiment. The data processing device 1 comprises a data acquisition unit 10, a controllability Gramian calculation unit 12 and an output unit 15. The data processing device 1 estimates the limit:









G

(

)




[

equation


140

]







of the controllability Gramian:









G

(
t
)




[

equation


138

]







defined by:










G

(
t
)

:=



0


t




e

A

τ



B


B




e


A



τ



d

τ






[

equation


2

]









in
:









t
=





[

equation


144

]







when:











x
˙

(
t
)

=


A


x

(
t
)


+

B


u

(
t
)







[

equation


1

]







holds, where:









x

(
t
)




[

equation


142

]







is an n-dimensional vector representing the state of the control object:









u

(
t
)




[

equation


139

]







is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix.


The data acquisition unit 10 acquires a set of time-series state data in a plurality of time intervals. In the following description, the data acquisition unit 10 acquires a set of time-series state data:










x

(


[


t
11

,

t

1

2



]

,

x

1

1



)

,

x

(


[


t

2

1


,

t

2

2



]

,

x

1

1



)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)





[

equation


5

]







for the following q time intervals.












[


t

i

1


,

t

i

2



]




(


i

=

1

,
2
,


,
q

)







[

equation


4

]







The controllability Gramian calculation unit 12 defines:










𝓏

(
t
)



R
n





[

equation


75

]







expressed as:









[

equation


74

]











𝓏
˙

(
t
)

=


A




𝓏

(
t
)






(
13
)







based on the sets of state data:









E

(
t
)




[

equation


145

]







acquired by the data acquisition unit 10, calculates:









[

equation


81

]













𝓏


(

t
2

)


X


𝓏

(

t
2

)


-



𝓏


(

t
1

)


X


𝓏

(

t
1

)



=

-




t
1




t
2






𝓏


(
t
)



BB




𝓏

(
t
)


dt







(
16
)













𝓏

(


t
+

t

i

1



,

t

i

1


,

x

i

1



)

=



(


E

(
t
)



E
0

-
1



)





x

i

1







[

equation


89

]







and estimates:










G

(

)

=
X




[

equation


146

]







by numerically obtaining the solution X of the following linear equation.














x

i

1



(


E

(
h
)



E
0

-
1



)




X

(


E

(
h
)



E
0

-
1



)





x

i

1



-


x

i

1





Xx

i

1




=

-



0


h




x

i

1



(


E

(
t
)



E
0

-
1



)










BB


(


E

(
t
)



E
0

-
1



)





x

i

1




dt

(


i
=
1

,
2
,


,
q

)






[

equation


6

]







The output unit 15 outputs the controllability Gramian estimated by the controllability Gramian calculation unit 12.


Here, the control objects are, for example, the following.


(Example 1) Power system. The time series data include the amount of electricity generated, the amount of electricity used, temperature and humidity, etc. at each power plant.


(Example 2) Robots. Time-series data include load, posture, external force and motion state, etc.


(Example 3) Human body. Time-series data include blood pressure, pulse rate, body temperature and blood composition, etc.


(Example 4) Chemical plant. Time-series data include temperature, humidity, atmospheric pressure and amount of dust in the air, etc.


(Example 5) Equipment to be maintained. Time-series data include failure alarm frequency, electrical resistance, metal fatigue and heat generation, etc.


These control objects are configured as network systems. The mathematical model describing the time-series data output from these systems are complex, and generally such a mathematical model itself is unknown. In contrast, by estimating and outputting controllability Gramians using this implementation, it is possible to determine which node of the control objects configured as network systems should be controlled as the input channel, while the mathematical model is still unknown. The above features are also common to all of the following forms.


According to this embodiment, the controllability Gramian can be estimated using a data-driven method based on a set of acquired continuous-driven state data for a system of which mathematical model is unknown.


The Second Embodiment

The set of time-series state data acquired by the data acquisition unit may include noise. In this case, the controllability Gramian calculation unit 12 estimates:










G

(

)

=
X




[

equation


146

]







by computing an approximate solution to Problem 1(i) from:









[

equation


108

]













x

i

1



(


E

(
h
)



E
0
+


)




X

(


E

(
h
)



E
0
+


)





x

i

1



-


x

i

1





Xx

i

1




=

-





0



h





x

i

1



(


E

(
t
)



E
0
+


)





BB


(


E

(
t
)



E
0
+


)





x

i

1




dt

(


i
=
1

,
2
,

,
p

)









(
23
)








instead of (20). Here:









E

(
t
)




[

equation


145

]









and
:









E
0




[

equation


86

]







are defined with respect to the p data and the following.









h

:=


min

i


{

1
,
2
,

,
p

}



(


t

i

2


-

t

i

1



)





[

equation


109

]







According to this embodiment, controllability Gramians can be estimated using a data-driven method even when the state data contains noise.


The Third Embodiment

When the state data contains noise, knowledge about the signs of some or all of the matrix components of the solution X of (23) may be given as prior knowledge. In this case, the controllability Gramian calculation unit 12 performs numerical calculations using such prior knowledge.


According to this method, controllability Gramians can be estimated more accurately using a data-driven method, even when the state data contains noise.


The Fourth Embodiment


FIG. 2 shows a functional block diagram of a data processing device 1 according to the fourth embodiment. The data processing device 2 comprises a data acquisition unit 10, a maximization condition calculation unit 14, and an output unit 15. The data processing device 2 estimates the matrix B which maximizes the trace:









tr

(

G

(

)

)




[

equation


141

]







of the limit:









G

(

)




[

equation


140

]







of a controllability Gramian:









G

(
t
)




[

equation


138

]







defined by:










G

(
t
)

:=





0



t




e

A

τ




BB




e


A



τ



d

τ






[

equation


2

]









in
:









t
=





[

equation


144

]







when:











x
.

(
t
)

=


Ax

(
t
)

+

Bu

(
t
)






[

equation


1

]







holds, where:









x

(
t
)




[

equation


142

]







is an n-dimensional vector representing the state of the control object:









u

(
t
)




[

equation


139

]







is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix.


The data acquisition unit 10 acquires sets of time-series state data in a plurality of time intervals. In the following description, the data acquisition unit 10 acquires a set of time-series state data:










x

(


[


t
11

,

t
12


]

,

x
11


)

,

x

(


[


t
21

,

t
22


]

,

x
11


)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)





[

equation


5

]







for the following q time intervals.












[


t

i

1



,

t

i

2



]




(


i
=
1

,
2
,


,
q

)







[

equation


4

]







The maximization condition calculation unit 14 estimates the matrix B which maximizes:









tr

(

G

(

)

)




[

equation


141

]







by numerically obtaining the solution:










Y
~

*




[

equation


9

]







of the following linear equation:













x


(


t

i

2


,

t

i

1


,

x

i

1



)



Yx

(


t

i

2


,

t

i

1


,

x

i

1



)


-


x

i

1





Yx

i

1




=

-






t

i

1






t

i

2







x


(

t
,

t

i

1


,

x

i

1



)



Qx

(

t
,

t

i

1


,

x

i

1



)



dt

(


i
=
1

,
2
,


,
q

)








[

equation


8

]







based on the sets of state data:










x

(


[


t
11

,

t
12


]

,

x
11


)

,

x

(


[


t
21

,

t
22


]

,

x
11


)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)





[

equation


5

]







acquired by the data acquisition unit 10.


The output unit 15 outputs the input matrix when the controllability Gramian is maximized based on the maximization condition estimated by the maximization condition calculation unit 14.


According to this embodiment, for a system of which mathematical model is unknown, the input matrix when the controllability Gramian is maximized can be estimated using a data-driven method based on a set of acquired continuous-time state data.


The Fifth Embodiment

In the fifth embodiment, the maximization condition calculation unit 14 obtains a first input matrix:










B
~

1
*




[

equation


11

]







by calculating the unit eigenvector corresponding to the maximum eigenvalue of:










Y
~

*




[

equation


9

]







when:









B







m
=
1








R

n
×
m







[

equation


47

]







holds.


According to this embodiment, effective input matrix can be obtained when:









B







m
=
1








R

n
×
m







[

equation


47

]







holds.


The Sixth Embodiment

In the sixth embodiment, the maximization condition calculation unit obtains a second input matrix:










B
~

2
*




[

equation


12

]







by calculating the n×n matrix in which the (k, k) component is 1 and the other components are 0, where the (k, k) component is the one that is the maximum among the diagonal components of:











Y
~

*

,




[

equation


9

]







when:









B
=

D
n





[

equation


49

]







holds.


According to this embodiment, effective input matrix can be obtained when:









B
=

D
n





[

equation


49

]







holds.


The Seventh Embodiment (1)


FIG. 3 shows a functional block diagram of the control system 3 according to the seventh embodiment. The control system 3 comprises a sensor 16, a data processing device 1 and a control unit 18. The data processing device 1 comprises a data acquisition unit 10 and a controllability Gramian calculation unit 12. Namely, the data processing device 1 is configured and operates the same as the data processing device 1 in FIG. 1. The control system 3 controls external control objects 100 by inputting control signals to the control objects 100.


The sensor 16 detects a set of time-series state data from the control objects 100. The sensor 16 then transmits the detected set of state data to the data acquisition unit 10 of the data processing device 1. The data processing device 1 estimates:










G

(

)

=
X




[

equation


146

]







and transmits it to the control unit 18. The control unit 18 generates control signals based on:










G

(

)

=
X




[

equation


146

]







and controls the control objects 100 using them.


According to this embodiment, it is possible to detect time-series state data from an external control object, to estimate the controllability Gramian using a data-driven method based on the detected set of state data and appropriately to control said control object based on the estimated controllability Gramian.


The Seventh Embodiment (2)


FIG. 19 is a functional block diagram of the control system 5 according to a variant of the seventh embodiment. Control system 5 comprises a sensor 16, data processing device 1 and a control unit 19. The data processing device 1 comprises a data acquisition unit 10 and a controllability Gramian calculation unit 12. The control unit 19 has a control input value determination unit 191 and an input channel determination unit 192. Namely, the control system 5 differs from the control system 3 in FIG. 3 in that the control unit 19 includes the control input value determination unit 191 and the input channel determination unit 192. Other configurations of the control system 5 are common to the control system 3. In the following, we will focus on the parts that differ from the control system 3 and omit redundant explanations.


The control input value determination unit 191 is sent the set of state data detected by the sensor 16. Based on this set of state data, the control input value determination unit 191 determines what control should be performed on the control object determined by the input channel determination unit 192. The control input value determination unit 191 sends the determined control contents to the input channel determination unit 192.


The input channel determination unit 192 is sent the estimated controllability Gramian. The input channel determination unit 192 determines which of the control objects 100 should be controlled based on this controllability Gramian. The input channel determination unit 192 selects the determined control object and executes the control determined by the control input value determination unit 191 on the selected control object.


According to this embodiment, it is possible to detect time-series state data from external control objects, to estimate controllability Gramians using a data-driven method based on a set of detected state data, to select an appropriate control object based on the estimated controllability Gramians and appropriately to control the selected control object.


The Eighth Embodiment


FIG. 4 is a functional block diagram of the control system 4 according to the eighth embodiment. The control system 4 comprises a sensor 16, a data processing device 2 and a control unit 18. The data processing device 2 comprises a data acquisition unit 10 and a maximization condition calculation unit 14. Namely, the data processing device 2 is configured the same as the data processing unit 2 in FIG. 2 and operates the same. The control system 4 controls external control objects 100 by inputting control signals to the control object 100.


The sensor 16 detects a set of time-series state data from the control objects 100. The sensor 16 transmits the detected set of state data to the data acquisition unit 10 of the data processing device 1. The data processing device 2 estimates the matrix B which maximizes:









tr

(

G

(

)

)




[

equation


141

]







and sends it to the control unit 18. The control unit 18 generates control signals based on the matrix B and uses them to control the control objects 100.


According to this embodiment, it is possible to detect time-series state data from external control objects, to estimate the input matrix when the controllability Gramian is maximum using a data-driven method based on the detected set of state data, and appropriately to control the control objects based on the estimated input matrix.


The Ninth Embodiment


FIG. 5 is a flowchart showing the processing procedure of the data processing method according to the ninth embodiment. This method estimates the limit:









G

(

)




[

equation


140

]







of the controllability Gramian:









G

(
t
)




[

equation


138

]







defined by:










G

(
t
)

:=



0
t



e

A

τ





BB




e

A


τ




d

τ







[

equation


2

]









in
:









t
=





[

equation


144

]







when:











x
˙

(
t
)

=


A


x

(
t
)


+

B


u

(
t
)







[

equation


1

]







holds, where:









x

(
t
)




[

equation


142

]







is an n-dimensional vector representing the state of the control object:









u

(
t
)




[

equation


139

]







is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix.


In step S1, the method obtains the time-series state data set:










x

(


[


t
11

,

t

1

2



]

,

x

1

1



)

,

x

(


[


t

2

1


,

t

2

2



]

,

x

1

1



)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)





[

equation


5

]







in the q time intervals:












[


t

i

1


,

t

i

2



]




(


i

=

1

,
2
,


,
q

)







[

equation


4

]







when:










u

(
t
)


0




[

equation


39

]







holds.


In step S2, the method defines:










𝓏

(
t
)





R
n





[

equation


75

]







expressed as:









(
13
)












𝓏
.

(
t
)

=


A




𝓏

(
t
)



,




[

equation


74

]







calculates:









(
16
)













𝓏


(

t
2

)


X


𝓏

(

t
2

)


-



𝓏


(

t
1

)


X


𝓏

(

t
1

)



=

-




t
1


t
2





𝓏



(
t
)



BB




𝓏

(
t
)


d

τ








[

equation


81

]













𝓏

(


t
+

t

i

1



,

t

i

1


,

x

i

1



)

=



(

E


(
t
)



E
0

-
1



)





x

i

1







[

equation


89

]







and estimates:










G

(

)

=
X




[

equation


146

]







with respect to:










E

(
t
)

:=




[

x

(


t
+

t
11


,

t
11

,

x
11


)





x


(


t
+

t
21


,

t
21

,

x
21


)









x


(


t
+

t

n

1



,

t

n

1


,

x

n

1



)


]








[

equation


157

]













E
0



:=
[




x

1

1





x

2

1








x

n

1





]





[

equation


158

]







by numerically obtaining the solution X of the following linear equation.












x

i

1



(


E

(
h
)



E
0

-
1



)




X

(

E


(
h
)



E
0

-
1



)





x

i

1



-


x

i

1





Xx

i

1







[

equation


6

]









=

-



0
h




x

i

1



(


E

(
t
)



E
0

-
1



)





BB


(


E

(
t
)



E
0

-
1



)





x

i

1



dt









(


1
=
1

,
2
,


,
q

)




In step S3, the method outputs the controllability Gramian estimated in step S2.


According to this embodiment, the controllability Gramian can be estimated by a computer using a data-driven method based on a set of acquired continuous-time state data for a system of which mathematical model is unknown.


The Tenth Embodiment


FIG. 6 is a flowchart showing the processing procedure of the data processing method according to the tenth embodiment.


This method estimates the matrix B which maximizes the trace:









t


r

(

G

(

)

)





[

equation


141

]







of the limit:









G

(

)




[

equation


140

]







of the controllability Gramian:









G

(
t
)




[

equation


138

]







defined by:










G

(
t
)


:
=





0



t




e

A

τ




BB




e



A





τ



d

τ






[

equation


2

]









in
:









t
=






[



equation







1



4



4




]








when:











x
˙

(
t
)

=


Ax

(
t
)

+

Bu


(
t
)







[

equation


1

]







holds, where:









x

(
t
)




[

equation


142

]







is an n-dimensional vector representing the state of the control object:










u
(
t

)




[

equation


139

]







is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix.


In step S1, the method obtains the time-series state data set:










x

(


[


t

1

1


,

t

1

2



]

,

x

1

1



)

,

x

(


[


t

2

1


,

t

2

2



]

,

x

1

1



)

,

,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)





[

equation


5

]







in the q time intervals:









[


t

i

1


,

t

i

2



]




[

equation


4

]









(


i

=

1

,
2
,

,
q

)




when:










u

(
t
)


0




[

equation


39

]







holds.


In step S4, the method estimates the matrix B which maximizes:









tr


(

G

(

)

)





[

equation


141

]







by numerically obtaining the solution:










Y
˜

*




[

equation


9

]







of the following linear equation:













x


(


t

i

2


,

t

i

1


,

x

i

1



)



Yx

(


t

i

2


,

t

i

1


,

x

i

1



)


-


x

i

1





Yx

i

1




=


-






t

i

1






t

i

2







x


(

t
,

t

i

1


,

x

i

1



)



Qx

(

t
,

t

i

1


,

x



i

1




)


dt







[

equation


8

]









(


i
=
1

,
2
,

,
q

)




when:









Q

:
=

I




[

equation


7

]







holds.


In step S5, the method outputs the input matrix when the controllability Gramian is maximized based on the maximization condition estimated in step S4.


According to this embodiment, the input matrix when the controllability Gramian is maximized can be estimated by a computer using a data-driven method based on a set of acquired continuous-time state data for a system of which mathematical model is unknown.


The Eleventh Embodiment

The eleventh embodiment is a program. This program estimates the limit:









G

(

)




[

equation


140

]







of the controllability Gramian:









G

(
t
)




[

equation


138

]







defined by:










G

(
t
)


:
=





0



t




e

A

τ




BB




e



A





τ



d

τ






[

equation


2

]







in:









t
=





[

equation


144

]







when:











x
˙

(
t
)

=


Ax

(
t
)

+

B


u

(
t
)







[

equation


1

]







holds, where:









x

(
t
)




[

equation


142

]







is an n-dimensional vector representing the state of the control object:









u

(
t
)




[

equation


139

]







is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix. This program make a computer perform


the step of obtaining the time-series state data set:










x

(


[


t
11

,

t

1

2



]

,


x

1

1



)

,

x

(


[


t
21

,

t
22


]

,

x

1

1



)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)





[

equation


5

]







in the q time intervals:












[


t

i

1


,

t

i

2



]




(


i

=

1

,
2
,


,
q

)







[

equation


4

]







when:










u

(
t
)


0




[

equation


39

]







holds,


the step of defining:










𝓏

(
t
)





R
n





[

equation


75

]







expressed as:









(
13
)












𝓏
.

(
t
)

=


A







𝓏

(
t
)




,




[

equation


74

]







calculating:









(
16
)













𝓏


(

t
2

)


X


𝓏

(

t
2

)


-



𝓏


(

t
1

)


X


𝓏

(

t
1

)



=

-




t
1


t
2





𝓏



(
t
)



BB




𝓏

(
t
)


dt







[

equation


81

]













𝒵

(


t
+

t

i

1



,

t

i

1


,


x

i

1



)

=



(


E

(
t
)



E
0

-
1



)





x

i

1







[

equation


89

]







and estimating:










G

(

)

=
X




[

equation


146

]







by numerically obtaining the solution X of the following linear equation:












x

i

1



(


E

(
h
)



E
0

-
1



)




X

(


E

(
h
)



E
0

-
1



)





x

i

1



-


x

i

1




X


x

i

1







[

equation


6

]









=

-



0
h




x

i

1



(


E

(
t
)



E
0

-
1



)


B




B


(


E

(
t
)



E
0

-
1



)





x

i

1



d

t









(


i
=
1

,
2
,



,
q

)




with respect to:











E

(
t
)

:=




[

x

(


t
+

t

1

1



,


t
11

,

x
11


)






x


(


t
+

t
21


,


t
21

,

x
21


)











x


(


t
+

t

n

1



,


t

n

1


,

x

n

1



)






]




[

equation


157

]














E
0

:

=

[




x
11




x

2

1








x

n

1





]





[

equation


158

]







and


the step of outputting the estimated controllability Gramian.


According to this embodiment, software can be implemented as a program to estimate the controllability Gramian using a data-driven method based on a set of acquired continuous-time state data for a system of which the mathematical model is unknown.


The Twelfth Embodiment

The twelfth embodiment is a program. This program estimates the matrix B which maximizes the trace:









t


r

(

G

(

)

)





[

equation


141

]







of the limit:









G

(

)




[

equation


140

]







of the controllability Gramian:









[

equation


138

]









G

(
t
)










defined by:









[

equation


2

]










G

(
t
)

:=



0


t




?


BB



?

d

τ















?

indicates text missing or illegible when filed




in:









[

equation


144

]









t
=











when:









[

equation


1

]











x
˙

(
t
)

=


Ax

(
t
)

+

Bu

(
t
)












holds, where:









[

equation


142

]









x

(
t
)










is an n-dimensional vector representing the state of the control object:









[

equation


139

]









u

(
t
)










is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix.


This program make a computer perform


the step of obtaining the time-series state data set:









[

equation


5

]










x

(


[


t
11

,

t
12


]

,

x
11


)

,

x

(


[


t
21

,

t
22


]

,

x
11


)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)











in the q time intervals:









[

equation


4

]










[


t

i

1


,

t

i

2



]




(


i

=

1

,
2
,


,
q

)











when:









[

equation


39

]










u

(
t
)


0










holds,


the step of estimating the matrix B which maximizes:









[

equation


141

]









tr

(

G

(

)

)










by numerically obtaining the solution:









[

equation


9

]










Y
˜

*










of the following linear equation:









[

equation


8

]






















x


(


t

i

2


,

t

i

1


,

x

i

1



)


Y


x

(


t

i

2


,

t

i

1


,

x

i

1



)


-


x

i

1




Y


x

i

1




=

-




t

i

1





t

i

2







x


(

t
,

t

i

1


,

x

i

1



)


𝒬


x

(

t
,

t

i

1


,

x

i

1



)


dt









(


i
=
1

,
2
,


,
q

)




when:









[

equation


7

]









𝒬
:=
I










holds and


the step of outputting the input matrix when the controllability Gramian is maximized based on the estimated maximization condition.


According to this embodiment, software can be implemented as a program to estimate the input matrix when the controllability Gramian is maximized by a computer using a data-driven method based on a set of acquired continuous-time state data for a system of which mathematical model is unknown.


[Verification 1]

The following shows the results by the first embodiment with respect to a set of noiseless time-series state data are presented. We obtain the solution of Problem 1(i) in the system shown in (1) when:









[

equation


147

]









n
=
3










holds. Here, we assume that the following hold.









[

equation


148

]









A
:=

[



0


1


0





-
1




-
2



5




0


0



-
1




]
















[

equation


149

]









B
:=

[



0




1




1



]











For comparison, the true values of controllability Gramians:










G

(

)

,




[

equation


140

]







which is obtained using (3) and “lyap” of Matlab™, are as follows.










G

(

)

=

[





6
.
8


1

3



0


0.875




0


2.438


0.875




0.875


0.875


0.5



]





[

equation


150

]








FIGS. 7 through 12 show the time-series data without noise. Specifically,



FIG. 7 shows the time-series data in the time interval:









[


t
11

,

t

1

2



]




[

equation


151

]








FIG. 8 shows the time-series data in the time interval:









[


t
21

,

t

2

2



]




[

equation


152

]








FIG. 9 shows the time-series data in the time interval:









[


t
31

,

t

3

2



]




[

equation


153

]








FIG. 10 shows the time-series data in the time interval:









[


t
41

,

t

4

2



]




[

equation


154

]








FIG. 11 shows the time-series data in the time interval:









[


t
51

,

t

5

2



]




[

equation


155

]








FIG. 12 shows the time-series data in the time interval:









[


t
61

,

t

6

2



]




[

equation


156

]







The data acquisition unit 10 of the data processing device 1 in FIG. 1 acquires time-series state data sets:










x

(


[


t

1

1


,

t

1

2



]

,


x

1

1



)

,

x

(


[


t

2

1


,

t

2

2



]

,


x

1

1



)

,

,

x

(


[


t

6

1


,

t

6

2



]

,


x

6

1



)





[

equation


159

]







in the six time intervals shown in FIGS. 7 through 12. In this case:









E

(
t
)




[

equation


145

]









and

:









E
0




[

equation


86

]







are determined by (17) and (18), respectively. Then, the following is obtained.










D
2

=

[




-
111.33



58.021



-
4.521



76.071


19.44


113.823





-
114.141




-
218.566




-
104.632




-
112.263




-
107.462



6.134




2.389


27.479



-
233.018



31.881



-
169.49



6.613





-
70.105




-
46.067




-
7.41




-
59.942




-
14.605



28.202





-
79.264




-
55.735




-
9.642




-
9.312



2.214


48.204





-
6.446



63.871



-
104.918



37.316



-
46.008



36.435



]





[

equation


160

]







This means the following.










rank



(

D
2

)


=

6
=


(

1
/
2

)



n

(

n
+
1

)







[

equation


161

]







Therefore, (22) holds and there exists a unique solution X to the linear equation (22).


The controllability Gramian calculation unit 12 solves (22) numerically to obtain the following values.









[

equation


162

]









X
=

[



6.813


0


0.875




0


2.438


0.875




0.875


0.875


0.5



]





(
31
)







This is consistent with the true value (30), indicating the usefulness of this embodiment.


[Verification 2]

The following shows the results by the second embodiment with respect to a set of time-series state data including noise. FIGS. 13 through 18 show the time-series data including noise. Specifically,



FIG. 13 shows the time-series data in the time interval:









[


t

1

1


,

t

1

2



]




[

equation


151

]








FIG. 14 shows the time-series data in the time interval:









[


t
21

,

t

2

2



]




[

equation


152

]








FIG. 15 shows the time-series data in the time interval:









[


t
31

,

t
32


]




[

equation


153

]








FIG. 16 shows the time-series data in the time interval:









[


t
41

,

t
42


]




[

equation


154

]








FIG. 17 shows the time-series data in the time interval:









[


t
51

,

t
52


]




[

equation


155

]








FIG. 18 shows the time-series data in the time interval:









[


t
61

,

t
62


]




[

equation


156

]







The signals in FIGS. 13 to 18 are the signals in FIGS. 7 to 12 plus noise with mean 0 and variance 0.1, respectively.


The controllability Gramian calculation unit 12 solves (23) numerically to obtain the following values.









[

equation


163

]









X
=

[



6.612


0.016


0.812




0.016


2.228


0.826




0.812


0.826


0.47



]





(
32
)







This result is largely consistent with the true value (30), although it is worse than (31).


[Verification 3]

The following is the actual result obtained using the third embodiment with respect to sets of time-series state data including noise. In this case, the sets of the time-series state data are also shown in FIGS. 13 through 18, as described above.


Here, we assume that we have prior knowledge that:









G

(

)




[

equation


140

]







has the following pattern of sign.









[



+


0


*




0


+


*




*


*


+



]




[

equation


164

]







Here, * indicates that the sign is unknown. Then, the problem of estimating:









G

(

)




[

equation


140

]







boils down to the following optimization problem.









[

equation


165

]









{





min

ξ


R
6









D
2


ξ

-

d
2




2









s
.
t
.


ξ
1


>
0

,


ξ
3

>
0

,


ξ
6

>
0

,


ξ
2

=
0









(
33
)







Here:









D
2



R

6
×
6






[

equation


166

]








and









d
2



R
6





[

equation


167

]







are related to (21) and used in (23).









ξ


R
6





[

equation


168

]







corresponds to:









rvec

(
X
)




[

equation


169

]







of (21).









ξ
i




[

equation


170

]







is the i-th component of:








ξ



[

equation


171

]












X
=

[



6.657


0


0.794




0


2.25


0.806




0.794


0.806


0.526



]





[

equation


172

]







is obtained by numerically computing (33). This approximate solution applying the above prior knowledge is closer to the true value (30) than the approximate solution (32).


[Verification 4]

The following describes the results obtained from the fourth embodiment using the data in FIGS. 7 through 12. In this case:










D
1

=

[




-
114.033



53.362



-
6.242



38.923



-
9.11




-
3.311






-
56.781




-
280.494




-
88.226




-
114.21




-
113.084




-
29.771





8.664


11.552



-
230.861



9.445



-
185.936




-
36.786






-
54.121




-
66.226




-
2.459




-
70.825




-
24.182




-
18.02






-
79.077




-
58.335




-
9.658




-
21.689




-
7.737




-
1.471






-
1.685



53.084



-
104.062



17.608



-
60.784




-
8.659




]





[

equation


173

]







is obtained. (12) holds since:











rank


{

D
1



)

=

6
=


(

1
/
2

)



n

(

n
+
1

)







[

equation


174

]







holds. Namely, there exists a unique solution Y for (10) with respect to:









Q
:=

I
.





[

equation


7

]







Numerical calculation in the embodiment yields the following.











Y
~

*

=

[





1
.
5


0

0





0
.
5


0

0





1
.
2


5

0







0
.
5


0

0





0
.
5


0

0





1
.
2


5

0







1
.
2


5

0



1.25




6
.
7


5

0




]





[

equation


175

]







[Verification 5]

The following is the results obtained from the fifth embodiment using the data in FIGS. 7 through 12. The first input matrix:











B
˜

1
*

=

[





0
.
2


2

3







0
.
1


9

2







0
.
9


5

6




]





[

equation


176

]







is obtained with respect to:









B
=




m
=
1




R

3
×
m







[

equation


48

]







by the numerical calculation according to the embodiment when:









B





m
=
1




R

n
×
m







[

equation


47

]







holds.









t


r

(

G

(

)

)





[

equation


141

]







takes the maximum value of:










t


r

(

G

(

)

)


=
7.293




[

equation


177

]







by this first input matrix. In contrast to this result:









t


r

(

G

(

)

)





[

equation


141

]







takes the value of:










tr


(

G

(

)

)


=
4.917




[

equation


179

]







when the input matrix:









B


R

3
×
3






[

equation


178

]







with all components ⅓ is used as an example for comparison. This is clearly smaller than the value:










t


r

(

G

(

)

)


=


7
.
2


9

3





[

equation


177

]







obtained in the embodiment.


[Verification 6]

The following are the results obtained from the sixth embodiment using the data in FIGS. 7 through 12. The second input matrix:











B
˜

2
*

=

[



0


0


0




0


0


0




0


0


1



]





[

equation


181

]







is obtained with respect to:










tr


(

G

(

)

)


=
6.75




[

equation


180

]







by the numerical calculation according to the embodiment when:









B
=

D
n





[

equation


49

]







holds.


The values of:









t


r

(

G

(

)

)





[

equation


141

]







in each case of:









B
=

D
n





[

equation


49

]







are shown in Table 1. From this, it is understood that:










B
~

2
*




[

equation


12

]







is the most effective input matrix, indicating the usefulness of this embodiment.












TABLE 1







B
tr(G(∞))



















[e1 0 0]
1.500



[0 e2 0]
0.500



{tilde over (B)}2* = [0 0 e3]
6.750










Embodiments in Case the Amount of Action Changes

In the embodiments described above, one of the goals was to estimate the “ease of control” index when the input channel to each node of the network is fixed (Problem 1(i)). In other words, the following problems were discussed.


(Problem 1)

For the system shown in Equation (1), we assume that the matrix A is a Hurwitz matrix and is unknown. We assume that a set of data:










x

(


[


t

i

1


,

t

i

2



]

,

x

i

1



)




(


i
=
1

,
2
,


,
q

)





[

equation


43

]







consisting of segments of q state trajectories is given.


Here:










t

i

1


<

t

i

2






[

equation


44

]







holds.


(i):









G

(

)




[

equation


140

]







is estimated when a known input matrix B:









x


(


[


t

i

1


,

t

i

2



]

,

x

i

1



)




(


i
=
1

,
2
,


,
q

)





[

equation


43

]







is given.


In other words, in the embodiments explained above, it was assumed that the amount of action between each node constituting the network, i.e., the strength of the connection between nodes and an indicator of the amount of transmitted data, are all fixed values. In this case, it was possible to know how to select input channels to improve ease of control. In contrast, the following describes the estimation of ease of control when the amount of action is adjusted, i.e., estimation of controllability Gramian when the amount of action changes. The goal is to know how to adjust the action amount in order to improve the ease of control.


In the following, we consider the estimation of the controllability Gramian when a change A is added to an unknown matrix A. This problem can be formulated as follows.


(Problem 3)

We consider the system in Equation (1). Here, we assume that:









B


R

n
×
m






[

equation


28

]







is known, while:









A


R

n
×
n






[

equation


27

]







is unknown. Also, we assume that the data of N state trajectories are given to the system as follows:









[

equation


182

]









:


x

(


[


t

i

1


,

t

i

2



]

,

x

i

1



)




(


i
=
1

,
2
,


,
N

)





(
34
)







Here:









t

i

1


<


t

i

2





(


i
=
1

,
2
,


,
N

)






[

equation


183

]







holds. Also, we assume that:









Δ


R

n
×
n






[

equation


185

]







such that:









A
+
Δ




[

equation


184

]







is a stable matrix is arbitrarily given. At this point, we obtain:










G
c

(

A
+
Δ

)




[

equation


186

]







Problem 3 is to estimate the value of the controllability Gramian from the data when the A matrix is varied by A for the system in Equation (1). This problem appears, for example, when applying state feedback to the system or when adjusting the connection strength between nodes in a network system.


[Data-Driven Solution of the Lyapunov Equation]

In the following, we describe a method for solving the Lyapunov equation using the state trajectories of the system as a preparation for Problem 3. We consider the following system:









[

equation


187

]











x
.

(
t
)

=


A
~



x

(
t
)






(
35
)







Here:






[

equation


188

]







A
~



R

n
×
n






is a stability matrix. Also, the state trajectory of the system in equation (35) can be expressed in the same way as described above as:






[

equation


189

]






x

(

t
,

t
1

,

x
1


)




For this system, we consider the following Lyapunov equation.






[

equation


190

]













(


A
~

+

Δ
~


)




Y

+

Y

(


A
~

+

Δ
~


)


=

-
Q






(
36
)








Here:






[

equation


191

]







Δ
~



R

n
×
n






is a matrix such that:






[

equation


192

]







A
~

+

Δ
~





is stable, and:






[

equation






193

]






Q


S
n





is a symmetric matrix. In this case, it is known that equation (36) has a unique solution.


At this point, we consider the following problem (Problem 4)


We consider the system in Equation (35). Here:






[

equation


194

]







A
~



R

n
×
n






is assumed to be stable but unknown. For this system, we assume that the data D of N state trajectories are given as in Problem 3. Also, we assume that:






[

equation


191

]








Δ
~




R

n
×
n






and
:







[

equation


193

]






Q


S
n





such that:






[

equation


192

]







A
~

+

Δ
~





is a stable matrix are arbitrarily given. In this case, we obtain the solution of Equation (36).


The solution to Problem 4 is given as follows. First, for both sides of equation (36), we consider the quadratic form with respect to the state trajectory:






[

equation






195

]






x

(
t
)




then:






[

equation


196

]













x


(
t
)



(




A
~




Y

+

Y


A
~



)



x

(
t
)


+



x


(
t
)



(




Δ
~




Y

+

Y


Δ
~



)



x

(
t
)



=


-


x


(
t
)




Qx

(
t
)






(
37
)







holds. Here:






[

equation


198

]












d
dt



(



x


(
t
)



Yx

(
t
)


)


+



x


(
t
)



(




Δ
~




Y

+

Y


Δ
~



)



x

(
t
)



=


-


x


(
t
)




Qx

(
t
)







(
38
)








is obtained since:






[

equation


197

]












x


(
t
)



(


Y


A
~


+



A
~




Y


)



x

(
t
)


=





x


(
t
)



Y

(


A
~



x

(
t
)


)


+



(


A
~



x

(
t
)


)





Yx

(
t
)









=





x


(
t
)


Y



x

.


(
t
)


+




x

.




(
t
)



Yx

(
t
)









=



d
dt



(



x


(
t
)



Yx

(
t
)


)









holds for the first term on the left side.


Furthermore:





[

equation


199

]













x


(

t
2

)


Y


x

(

t
2

)


-



x


(

t
1

)


Y


x

(

t
1

)


+




t

1




t
2





x


(
t
)



(




Δ
~




Y

+

Y


Δ
~



)



x

(
t
)


dt



=


-




t
1


t
2





x


(
t
)



Qx

(
t
)


dt








(
39
)








is obtained by integrating both sides of this equation over the time interval:






[

equation


41

]







[


t
1

,

t
2


]

.




Here, N linear equations for the unknown number Y:









[

equation


201

]











(
40
)














x
i


(

t

i

2


)


Y



x
i

(

t

i

2


)


-



x
i


(

t

i

1


)


Y



x
i

(

t

i

1


)


+




t

i

1





t

i

2







x
i


(
t
)



(




Δ
~




Y

+

Y


Δ
~



)




x
i

(
t
)



dt



=

-




t

i

1





t

i

2







x
i


(
t
)


𝒬



x
i

(
t
)



dt









(


i
=
1

,
2
,


,
N

)




is obtained by applying the state trajectory data:









[

equation


200

]











x
i

(
t
)

:=

x

(

t
,

t

i

1


,

x

i

1



)











to equation (39).


Next, we consider rewriting equation (40) in a matrix-based form. First, equation (39) can be expressed as:









[

equation


202

]











(
41
)












(


lvec

(


x

(

t
2

)

,

x

(

t
2

)


)

-

lvec

(


x

(

t
1

)

,

x

(

t
1

)


)

+




t
1




t
2





(


lvec

(



Δ
~



x

(
t
)


,

x

(
t
)


)

+

lvec

(


x

(
t
)

,


Δ
~



x

(
t
)



)


)


dt



)




rvec

(
Y
)


=

-




t
1




t
2





lvec

(


x

(
t
)

,

x

(
t
)


)


dt



rvec

(
𝒬
)








using lvec and rvec. If this is noted, Equation (40) can be rewritten as:









[

equation


203

]












D
~

1




rvec

(
Y
)


=



D
~

2




rvec

(
𝒬
)






(
42
)







equivalently. Here:









[

equation


204

]











D
~

1

,



D
~

2



R

N
×

(

1
/
2

)



n

(

n
+
1

)














are matrices determined by the state trajectory data of the system in Equation (35) and their i-th row:









[

equation


205

]











d
~


1

i


,



d
~


2

i




R

1
×

(

1
/
2

)



n

(

n
+
1

)














are given by:









[

equation


206

]



















?

:=


lvec

(



x
i

(

t

i

2


)

,

x

(

t

i

2


)


)

-

lvec

(



x
i

(

t

i

1


)

,


x
i

(

t

i

1


)


)

+




t

i

1





t

i

2






(


levc

(



?



x
i

(
t
)


,


x
i

(
t
)


)

+

lvec

(



x
i

(
t
)

,


Δ
~




x
i

(
t
)



)


)



dt












[

equation


207

]











d
~


2

i


:=

-




t

i

1





t

i

2






lvec

(



x
i

(
t
)

,


x
i

(
t
)


)



dt
















?

indicates text missing or illegible when filed




respectively.


The solution to Problem 4 is then obtained as follows.


(Theorem 3)

We consider Problem 4. If:









[

equation


68

]










rank
(

D
1

)

=


1
2



n

(

n
+
1

)






(
12
)







holds, the solution of equation (40) is unique in the range of symmetric matrices and is equal to the solution of Problem 4. (End of Theorem 3).


[Data-Driven Estimation of Controllability Gramians].

Next, we consider the solution to Problem 3.


Data Transformation








[

equation


186

]










G
c

(

A
+
Δ

)










can be computed using Theorem 3 if state trajectory data of:









[

equation


74

]











z
˙

(
t
)

=


A




z

(
t
)






(
13
)







which is a dual system of Equation (35) are available. This proposition can be expressed as follows. From Lemma 1, the controllability:









[

equation


186

]










G
c



(

A
+
Δ

)











is obtained as a unique solution of where equation of:









[

equation


208

]












(

A
+
Δ

)


X

+


X

(

A
+
Δ

)




=

-


BB


.







(
43
)








Similarly to the above discussion:









[

equation


212

]











(
44
)














z


(

t
2

)


X


z

(

t
2

)


-



z


(

t
1

)


X


z

(

t
1

)


+


?



z


(
t
)



(


Δ

X

+

X


Δ




)




z

(
t
)



dt


=

-




t
1




t
2






z


(
t
)


BB


z

(
t
)



dt










?

indicates text missing or illegible when filed




is obtained if we note that equation (43) corresponds to equation (36) to which:









[

equation


209

]










A
~

=

A

















[

equation


210

]










Δ
~

=

Δ

















[

equation


211

]









𝒬
=

BB












are applied. Thus, by applying the state trajectory of the system in equation (13) to equation (44), a linear equation with the solution of:









[

equation


186

]










G
c

(

A
+
Δ

)










can be constructed from the state trajectory data only.


However, in Problem 3, the state trajectory of the system in Equation (13) is not generally available except in the case that









A
=

A






[

equation


213

]







holds. Therefore, the above ideas cannot be used directly for Problem 3. On the other hand, if the state trajectory of the system in Equation (35) can be transformed into the state trajectory of the system in Equation (13), a solution to Problem 3 can be constructed. Such a data transformation is accomplished by the following Lemma 5.


(Lemma 5)

We express the state trajectory of the system of Equation (13) by:









z

(

t
,

t

i

1


,

z

i

1



)




[

equation


214

]







in the same way as aforementioned









Q



S
n

.





[

equation


193

]







We assume that the state trajectory data D in Problem 3 is arbitrarily given. Here, we assume that N≥n. We assume that, for this D:










E

(
t
)



R

n
×
n






[

equation


215

]













E
0



R

n
×
n






[

equation


216

]












h


(

0
,


)





[

equation


217

]







are defined by:









[

equation


218

]










E

(
t
)

:=

[


x

(


t
+

t
11


,

t
11

,

x
11


)



x

(


t
+

t
21


,

t
21

,

x
21


)







x

(


t
+

t

n

1



,

t

n

1


,

x

n

1



)


]






(
45
)











E
0

:=

[




x
11




x
21







x

n

1





]







h
:=


min

i


{

1
,
2
,

,
n

}



(


t

i

2


-

t

i

1



)





respectively. In this case, if E is regular:









[

equation


221

]










z

(

t
,

t

i

1


,

x

i

1



)

=



(


E

(
t
)



E
0

-
1



)





x

i

1







(
46
)







holds for arbitrary:









i


{

1
,
2
,


,
N

}





[

equation


219

]









and
:









t



[

0
,
h

]

.





[

equation


220

]







[Data-Driven Estimation].

We give the solution to problem 3 using the above results. First, we assume that the function:











z
i

(
t
)




R
n

(


i
=
1

,
2
,


,
N

)





[

equation


222

]







is defined by the following.











z
i

(
t
)

:=



(


E

(
t
)



E
0

-
1



)





x

i

1







[

equation


223

]







Then N linear equations:









[

equation


225

]













z
i


(
h
)




Xz
i

(
h
)


-



z
i


(
0
)




Xz
i

(
0
)


+



0
h




z
i


(
t
)



(


Δ

X

+

X


Δ




)




z
i

(
t
)


dt



=

-



0
h




z
i


(
t
)




BBz
i

(
t
)



dt

(


i
=
1

,
2
,


,
N

)








(
47
)







are obtained by applying:










z

(
t
)

=



z
i

(
t
)



(


i
=
1

,
2
,


,
N

)






[

equation


224

]







to equation (44). Furthermore:











D
1



rvec

(
X
)


=


D
2



rvec

(

BB


)






[

equation


226

]







is obtained as an equivalent expression to equation (47). is obtained. Here:










D
1

,


D
2



R

N
×

(

1
/
2

)



n

(

n
+
1

)








[

equation


227

]







are matrices determined by the state trajectory data D in Problem 3 and their i-th row:










d

1

i


,


d

2

i




R

1
×

(

1
/
2

)



n

(

n
+
1

)








[

equation


228

]







are given by:










d

1

i


:=


lvec

(



z
i

(
h
)

,


z
i

(
h
)


)

-

lvec

(



z
i

(
0
)

,


z
i

(
0
)


)

+



0
h



(


lvec

(



Δ





z
i

(
t
)


,


z
i

(
t
)


)

+

lvec

(



z
i

(
t
)

,


Δ





z
i

(
t
)



)


)


dt







[

equation


229

]













d

2

i


:=

-



0
h



lvec

(



z
i

(
t
)

,


z
i

(
t
)


)


dt







[

equation


230

]







respectively.


Thus, the solution to Problem 3 can be obtained as follows.


(Theorem 4)

We consider Problem 3. Here, we assume that N≥n. If E0 in Equation (45) is regular and:









rank



(

D
1

)

=


1
2



n

(

n
+
1

)







[

equation


68

]







holds, the solution of Equation (47) is unique in the range of symmetric matrices and is equal to the solution of Problem 3.


The Thirteenth Embodiment

As in the first embodiment, FIG. 1 is used to describe the data processing device according to the thirteenth embodiment.



FIG. 1 is a functional block diagram of the data processing device 1 according to the thirteenth embodiment. The data processing device 1 comprises a data acquisition unit 10, a controllability Gramian calculation unit 12 and an output unit 15.


The data processing device 1 estimates:










G
c

(

A
+
Δ

)




[

equation


232

]







when:











x
˙

(
t
)

=


A


x

(
t
)


+

B


u

(
t
)







[

equation


1

]







holds and the limit of controllability Gramian of matrix A, which is defined by:











G

(
t
)

:=



0
t



c

A

τ




BB




c


A



τ



d

τ



,




[

equation


2

]







is defined by:










G
c

(
A
)




[

equation


231

]







where:









x

(
t
)




[

equation


142

]







is an n-dimensional vector representing the state of the control object:









u

(
t
)




[

equation


139

]







is an m-dimensional vector representing the control input, A is an unknown n×n matrix and B is a known n×m matrix.


The data acquisition unit 10 acquires a set of time-series state data in a plurality of time intervals. In the following description, the data acquisition unit 10 acquires a set of time-series state data:










x

(


[


t
11

,

t

1

2



]

,

x

1

1



)

,

x

(


[


t

2

1


,

t

2

2



]

,

x

1

1



)

,


,

x

(


[


t

q

1


,

t

q

2



]

,

x

q

1



)





[

equation


5

]







for the following q time intervals.










[


t

i

1


,

t

i

2



]




(


i

=

1

,
2
,


,
q

)





[

equation


4

]







The controllability Gramian calculation unit 12 defines:










z

(
t
)



R
n





[

equation


75

]







expressed as:









[

equation


74

]











z
.

(
t
)

=


A




z

(
t
)







(
13
)








based on the sets of state data:









E

(
t
)




[

equation


145

]







acquired by the data acquisition unit 10, calculates:









[

equation


81

]













z


(

t
2

)



Xz

(

t
2

)


-



z


(

t
1

)



Xz

(

t
1

)



=

-




t
1


t
2





z


(
t
)



BB




z

(
t
)


dt







(
16
)













z

(


t
+

t

i

1



,

t

i

1


,

x

i

1



)

=



(


E

(
t
)



E
0

-
1



)





x

i

1







[

equation


89

]







and estimates:











G
c

(

A
+
Δ

)

=
X




[

equation


233

]







by numerically obtaining the solution X of the following linear equation.









[

equation


225

]













z
i


(
h
)




Xz
i

(
h
)


-



z
i


(
0
)




Xz
i

(
0
)


+



0
h




z
i


(
t
)



(


Δ

X

+

X


Δ




)




z
i

(
t
)


dt



=

-



0
h




z
i


(
t
)


BB



z
i

(
t
)


dt



(


i
=
1

,
2
,


,
N

)








(
47
)







The output unit 15 outputs the controllability Gramian estimated by the controllability Gramian calculation unit 12.


According to this embodiment, the controllability Gramian can be estimated when the amount of action between each node constituting the network changes.


The present disclosure has been described above based on the embodiments. It is understood by those skilled in the art that these embodiments are examples, that various variations are possible in the combination of each component and each processing process and that such variations are also within the scope of the disclosure.


Any combination of the above mentioned embodiments and variations is also useful as an embodiment of the disclosure. The new embodiment resulting from the combination will have the respective effects of each of the embodiments and variations that are combined.


INDUSTRIAL APPLICABILITY

The principle of the disclosure can be applied to the control of systems in various fields as follows.


(Application example 1) When the control object is an electric power system, optimal power supply control can be performed based on time-series data of power generation, power consumption, temperature and humidity, etc. at each power plant.


(Application example 2) When the control object is a robot, optimal posture control can be performed based on time-series data such as load, posture, external force and motion state, etc.


(Application example 3) When the control object is a human body, optimal medication control can be performed based on time-series data such as blood pressure, pulse rate, body temperature and blood composition, etc.


(Application example 4) When the control object is a chemical plant, the system can perform optimal operation control based on time-series data such as temperature, humidity, air pressure and the amount of dust in the air, etc.


(Application example 5) When the control object is equipment maintenance, optimal maintenance control can be performed based on time-series data such as failure alarm frequency, electrical resistance, metal fatigue and heat generation, etc.


REFERENCE SIGNS LIST






    • 1. Data processing device.


    • 2. Data processing device.


    • 3. Control system.


    • 4. Control system.


    • 5. Control system.


    • 10. Data acquisition unit.


    • 12. Controllability Gramian calculation unit.


    • 14. Maximization condition calculation unit.


    • 15. Output unit.


    • 16. Sensor.


    • 18. Control unit.


    • 19: Control unit.


    • 100: Control object.


    • 191: Control input value determination unit.


    • 192: Input channel determination unit.

    • S1: The step of acquiring a set of time-series state data.

    • S2: The step of estimating the controllability Gramian.

    • S3: The step of outputting the controllability Gramian.

    • S4: The step of estimating the condition when the controllability Gramian is maximum.

    • S5: The step of outputting the input matrix when the controllability Gramian is maximized.




Claims
  • 1-13. (canceled)
  • 14. A data processing device to estimate the limit:
  • 15. The data processing device according to claim 14, wherein the set of time-series state data acquired by the data acquisition unit includes noise andwherein the controllability Gramian calculation unit estimates:
  • 16. The data processing device according to claim 15, wherein the controllability Gramian calculation unit performs numerical calculations using prior knowledge about the signs of some or all of the matrix components of the solution X.
  • 17. A data processing device to estimate the matrix B which maximizes the trace:
  • 18. The data processing device according to claim 17, wherein the maximization condition calculation unit obtains a first input matrix:
  • 19. The data processing device according to claim 17, wherein the maximization condition calculation unit obtains a second input matrix:
  • 20. A control system to control external control objects, comprising: a sensor that detects a set of time-series state data from the control objects;
  • 21. A control system to control external control objects, comprising: a sensor that detects a set of time-series state data from the control objects;the data processing device according to claim 14; anda control unit to control the control objects,wherein the sensor transmits the detected set of state data to the data acquisition unit of the data processing device,wherein the data processing device transmits the estimated matrix B which maximizes:
  • 22. A data processing method to estimate the limit:
  • 23. A data processing method to estimate the matrix B which maximizes the trace:
  • 24. A non-transitory computer readable medium that stores a program to estimate the limit:
  • 25. A non-transitory computer readable medium that stores a program to estimate the matrix B which maximizes the trace:
  • 26. A data processing device to estimate:
Priority Claims (1)
Number Date Country Kind
2021-143668 Sep 2021 JP national
PCT Information
Filing Document Filing Date Country Kind
PCT/JP2022/033059 9/2/2022 WO