In [1]:
import numpy as np
import control as ct
from numpy.linalg import eig
from numpy import linalg as LA
In [2]:
A = np.asarray([
    [0, 1, 0],
    [0, 0, 1],
    [-6, -11, -6]
])
B = np.asarray([
    [0],
    [0],
    [6]
])
C = np.asarray([1, 0, 0])
D = np.asarray([0])

sys = ct.ss(A, B, C, D)
In [3]:
sys
Out[3]:
$$ \left(\begin{array}{rllrllrll|rll} 0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ 0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ -6\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&-11\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&-6\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&6\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ \hline 1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ \end{array}\right) $$
In [4]:
ct.ss2tf(sys)
Out[4]:
$$\frac{-1.776 \times 10^{-15} s^2 - 5.329 \times 10^{-15} s + 6}{s^3 + 6 s^2 + 11 s + 6}$$
In [5]:
eigenvalues, eigenvectors = LA.eig(A)
inv_eigen = LA.inv(eigenvectors)

Am = (inv_eigen @ A) @ eigenvectors
Bm = inv_eigen @ B
Cm = C @ eigenvectors
In [6]:
sys2 = ct.ss(Am, Bm, Cm, D)
sys2
Out[6]:
$$ \left(\begin{array}{rllrllrll|rll} -1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&2.&\hspace{-1em}44&\hspace{-1em}\cdot10^{-15}&-4.&\hspace{-1em}11&\hspace{-1em}\cdot10^{-15}&-5.&\hspace{-1em}2&\hspace{-1em}\phantom{\cdot}\\ -3.&\hspace{-1em}55&\hspace{-1em}\cdot10^{-15}&-2\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&-7.&\hspace{-1em}11&\hspace{-1em}\cdot10^{-15}&-27.&\hspace{-1em}5&\hspace{-1em}\phantom{\cdot}\\ -5.&\hspace{-1em}33&\hspace{-1em}\cdot10^{-15}&3.&\hspace{-1em}55&\hspace{-1em}\cdot10^{-15}&-3\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&-28.&\hspace{-1em}6&\hspace{-1em}\phantom{\cdot}\\ \hline -0.&\hspace{-1em}577&\hspace{-1em}\phantom{\cdot}&0.&\hspace{-1em}218&\hspace{-1em}\phantom{\cdot}&-0.&\hspace{-1em}105&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\ \end{array}\right) $$
In [7]:
ct.ss2tf(sys2)
Out[7]:
$$\frac{-1.776 \times 10^{-15} s^2 - 7.105 \times 10^{-15} s + 6}{s^3 + 6 s^2 + 11 s + 6}$$
In [8]:
eigenvectors
Out[8]:
array([[-0.57735027,  0.21821789, -0.10482848],
       [ 0.57735027, -0.43643578,  0.31448545],
       [-0.57735027,  0.87287156, -0.94345635]])

Example 9-5¶

https://dynamics-and-control.readthedocs.io/en/latest/0_Getting_Started/Cheatsheet.html#id2

In [9]:
!pip install sympy

Requirement already satisfied: sympy in /usr/local/lib/python3.9/dist-packages (1.12)
Requirement already satisfied: mpmath>=0.19 in /usr/local/lib/python3.9/dist-packages (from sympy) (1.3.0)
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv
In [10]:
import sympy
sympy.init_printing()
s = sympy.Symbol('s')  # A single symbol
In [11]:
A = np.asarray([
    [0,   1],
    [-2, -3]])
identity_matrix = np.identity(2)
laplce_target = sympy.Matrix(s * identity_matrix -A)
In [12]:
laplce_target_inv = laplce_target.inv()
In [13]:
laplce_target_inv
Out[13]:
$\displaystyle \left[\begin{matrix}\frac{2.0 s + 6.0}{2.0 s^{2} + 6.0 s + 4.0} & - \frac{2}{- 2.0 s^{2} - 6.0 s - 4}\\\frac{2}{- 1.0 s^{2} - 3.0 s - 2} & \frac{1.0 s}{1.0 s^{2} + 3.0 s + 2.0}\end{matrix}\right]$
In [14]:
type(laplce_target_inv)
Out[14]:
sympy.matrices.dense.MutableDenseMatrix
In [15]:
t, s = sympy.symbols('t, s')
a, b, c, d = sympy.symbols('a b c d')
partial = sympy.matrices.dense.MutableDenseMatrix(2,2, [a, b, c, d])

values = {a: laplce_target_inv[0].apart(s),
          b: laplce_target_inv[1].apart(s),
          c: laplce_target_inv[2].apart(s),
          d: laplce_target_inv[3].apart(s)
         }
partial = partial.subs(values)
In [16]:
partial
Out[16]:
$\displaystyle \left[\begin{matrix}\frac{2.0}{s + 1} - \frac{0.5}{0.5 s + 1.0} & \frac{1.0}{s + 1} - \frac{0.5}{0.5 s + 1.0}\\- \frac{2.0}{s + 1} + \frac{1.0}{0.5 s + 1.0} & - \frac{1.0}{s + 1} + \frac{1.0}{0.5 s + 1.0}\end{matrix}\right]$
In [17]:
type(partial)
Out[17]:
sympy.matrices.dense.MutableDenseMatrix
In [18]:
a, b, c, d = sympy.symbols('a b c d')
complete = sympy.matrices.dense.MutableDenseMatrix(2,2, [a, b, c, d])

values = {a: sympy.inverse_laplace_transform(partial[0], s, t),
          b: sympy.inverse_laplace_transform(partial[1], s, t),
          c: sympy.inverse_laplace_transform(partial[2], s, t),
          d: sympy.inverse_laplace_transform(partial[3], s, t)
         }

complete = complete.subs(values)
In [19]:
complete
Out[19]:
$\displaystyle \left[\begin{matrix}- 1.0 e^{- 2.0 t} \theta\left(t\right) + 2.0 e^{- t} \theta\left(t\right) & - 1.0 e^{- 2.0 t} \theta\left(t\right) + 1.0 e^{- t} \theta\left(t\right)\\2.0 e^{- 2.0 t} \theta\left(t\right) - 2.0 e^{- t} \theta\left(t\right) & 2.0 e^{- 2.0 t} \theta\left(t\right) - 1.0 e^{- t} \theta\left(t\right)\end{matrix}\right]$

3.3. What is that θ?

The unit step function is also known as the Heaviside step function. We will see this function often in inverse laplace transforms. It is typeset as by sympy.

In [20]:
sympy.Heaviside(t)
Out[20]:
$\displaystyle \theta\left(t\right)$
In [21]:
#예제
a = sympy.symbols('a', real=True, positive=True)
t, s = sympy.symbols('t, s')
f = sympy.exp(-a*t)
sympy.laplace_transform(f, t, s)
Out[21]:
$\displaystyle \left( \frac{1}{a + s}, \ - a, \ \text{True}\right)$

A matrix를 Diagonal form으로 변경했을 경우 계산¶

In [22]:
A
eigenvalues, eigenvectors = LA.eig(A)
inv_eigen = LA.inv(eigenvectors)
t, s = sympy.symbols('t, s')

print(eigenvalues)
print(eigenvectors)
diagonal_matrix = inv_eigen @ A @ eigenvectors
lamda_matrix = sympy.diag(sympy.exp(eigenvalues[0]*t),sympy.exp(eigenvalues[1]*t) )
# P^-1 @ A P가 아니라 P @ A @ P^-1인 형식
complete2 = eigenvectors @ lamda_matrix @ inv_eigen

[-1. -2.]
[[ 0.70710678 -0.4472136 ]
 [-0.70710678  0.89442719]]
In [23]:
lamda_matrix
Out[23]:
$\displaystyle \left[\begin{matrix}e^{- 1.0 t} & 0\\0 & e^{- 2.0 t}\end{matrix}\right]$
In [24]:
complete2
Out[24]:
$\displaystyle \left[\begin{matrix}- 1.0 e^{- 2.0 t} + 2.0 e^{- 1.0 t} & - 1.0 e^{- 2.0 t} + 1.0 e^{- 1.0 t}\\2.0 e^{- 2.0 t} - 2.0 e^{- 1.0 t} & 2.0 e^{- 2.0 t} - 1.0 e^{- 1.0 t}\end{matrix}\right]$

Example 9-6¶

In [25]:
A = np.asarray([
    [0,  1],
    [-2, -3]
])
B = np.asarray([
    [0],
    [1]
])
s = sympy.symbols('s', real=True, positive=True)
identity_matrix = np.identity(2)
Pi_matrix = sympy.Matrix(s * identity_matrix - A)
Pi_matrix = Pi_matrix.inv()
initial_tmp = Pi_matrix
response_tmp = Pi_matrix @ B
In [26]:
a, b, c, d = sympy.symbols('a b c d')
initial = sympy.matrices.dense.MutableDenseMatrix(2,2, [a, b, c, d])
values = {a: Pi_matrix[0].apart(s),
          b: Pi_matrix[1].apart(s),
          c: Pi_matrix[2].apart(s),
          d: Pi_matrix[3].apart(s)
         }
initial = initial.subs(values)

response = sympy.matrices.dense.MutableDenseMatrix(2,1, [a, b])
values = {a: response_tmp[0].apart(s),
          b: response_tmp[1].apart(s),
         }
response = response.subs(values)
In [27]:
initial
Out[27]:
$\displaystyle \left[\begin{matrix}\frac{2.0}{s + 1} - \frac{0.5}{0.5 s + 1.0} & \frac{1.0}{s + 1} - \frac{0.5}{0.5 s + 1.0}\\- \frac{2.0}{s + 1} + \frac{1.0}{0.5 s + 1.0} & - \frac{1.0}{s + 1} + \frac{1.0}{0.5 s + 1.0}\end{matrix}\right]$
In [28]:
response
Out[28]:
$\displaystyle \left[\begin{matrix}\frac{1.0}{s + 1} - \frac{0.5}{0.5 s + 1.0}\\- \frac{1.0}{s + 1} + \frac{1.0}{0.5 s + 1.0}\end{matrix}\right]$
In [29]:
a, b, c, d = sympy.symbols('a b c d')
#t, s = sympy.symbols('t, s')

initial_time = sympy.matrices.dense.MutableDenseMatrix(2,2, [a, b, c, d])

values = {
    a: sympy.inverse_laplace_transform(initial[0], s, t),
    b: sympy.inverse_laplace_transform(initial[1], s, t),
    c: sympy.inverse_laplace_transform(initial[2], s, t),
    d: sympy.inverse_laplace_transform(initial[3], s, t)
}
initial_time = initial_time.subs(values)

response_time = sympy.matrices.dense.MutableDenseMatrix(2,1, [a, b])
values = {
    a: sympy.inverse_laplace_transform(response[0], s, t),
    b: sympy.inverse_laplace_transform(response[1], s, t)

}
response_time = response_time.subs(values)
In [30]:
initial_time
Out[30]:
$\displaystyle \left[\begin{matrix}- 1.0 e^{- 2.0 t} \theta\left(t\right) + 2.0 e^{- t} \theta\left(t\right) & - 1.0 e^{- 2.0 t} \theta\left(t\right) + 1.0 e^{- t} \theta\left(t\right)\\2.0 e^{- 2.0 t} \theta\left(t\right) - 2.0 e^{- t} \theta\left(t\right) & 2.0 e^{- 2.0 t} \theta\left(t\right) - 1.0 e^{- t} \theta\left(t\right)\end{matrix}\right]$
In [31]:
response_time
Out[31]:
$\displaystyle \left[\begin{matrix}- 1.0 e^{- 2.0 t} \theta\left(t\right) + 1.0 e^{- t} \theta\left(t\right)\\2.0 e^{- 2.0 t} \theta\left(t\right) - 1.0 e^{- t} \theta\left(t\right)\end{matrix}\right]$

response_time에서 t를 t-tau 로 변경하고, 적분