8000 Add similarity transformation function (#319) · Gonewithwind0/python-control@fbefeba · GitHub
[go: up one dir, main page]
Skip to content

Commit fbefeba

Browse files
murrayrmmurrayrm
authored
Add similarity transformation function (python-control#319)
* add similarity transformation function (w/ optional time rescaling) * update unit test to avoid bugs in modal form
1 parent 732149c commit fbefeba

File tree

2 files changed

+109
-4
lines changed

2 files changed

+109
-4
lines changed

‎control/canonical.py

Lines changed: 37 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,11 @@
66
from .statesp import StateSpace
77
from .statefbk import ctrb, obsv
88

9-
from numpy import zeros, shape, poly, iscomplex, hstack
9+
from numpy import zeros, shape, poly, iscomplex, hstack, dot, transpose
1010
from numpy.linalg import solve, matrix_rank, eig
1111

12-
__all__ = ['canonical_form', 'reachable_form', 'observable_form']
12+
__all__ = ['canonical_form', 'reachable_form', 'observable_form', 'modal_form',
13+
'similarity_tra 10000 nsform']
1314

1415
def canonical_form(xsys, form='reachable'):
1516
"""Convert a system into canonical form
@@ -212,3 +213,37 @@ def modal_form(xsys):
212213
zsys.C = xsys.C.dot(Tzx)
213214

214215
return zsys, Tzx
216+
217+
218+
def similarity_transform(xsys, T, timescale=1):
219+
"""Perform a similarity transformation, with option time rescaling.
220+
221+
Transform a linear state space system to a new state space representation
222+
z = T x, where T is an invertible matrix.
223+
224+
Parameters
225+
----------
226+
T : 2D invertible array
227+
The matrix `T` defines the new set of coordinates z = T x.
228+
timescale : float
229+
If present, also rescale the time unit to tau = timescale * t
230+
231+
Returns
232+
-------
233+
zsys : StateSpace object
234+
System in transformed coordinates, with state 'z'
235+
236+
"""
237+
# Create a new system, starting with a copy of the old one
238+
zsys = StateSpace(xsys)
239+
240+
# Define a function to compute the right inverse (solve x M = y)
241+
def rsolve(M, y):
242+
return transpose(solve(transpose(M), transpose(y)))
243+
244+
# Update the system matrices
245+
zsys.A = rsolve(T, dot(T, zsys.A)) / timescale
246+
zsys.B = dot(T, zsys.B) / timescale
247+
zsys.C = rsolve(T, zsys.C)
248+
249+
return zsys

‎control/tests/canonical_test.py

Lines changed: 72 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,9 @@
22

33
import unittest
44
import numpy as np
5-
from control import ss, tf
5+
from control import ss, tf, tf2ss, ss2tf
66
from control.canonical import canonical_form, reachable_form, \
7-
observable_form, modal_form
7+
observable_form, modal_form, similarity_transform
88
from control.exception import ControlNotImplemented
99

1010
class TestCanonical(unittest.TestCase):
@@ -218,6 +218,76 @@ def test_arguments(self):
218218
np.testing.assert_raises(
219219
ControlNotImplemented, canonical_form, sys, 'unknown')
220220

221+
def test_similarity(self):
222+
"""Test similarty transform"""
223+
224+
# Single input, single output systems
225+
siso_ini = tf2ss(tf([1, 1], [1, 1, 1]))
226+
for form in 'reachable', 'observable':
227+
# Convert the system to one of the canonical forms
228+
siso_can, T_can = canonical_form(siso_ini, form)
229+
230+
# Use a similarity transformation to transform it back
231+
siso_sim = similarity_transform(siso_can, np.linalg.inv(T_can))
232+
233+
# Make sure everything goes back to the original form
234+
np.testing.assert_array_almost_equal(siso_sim.A, siso_ini.A)
235+
np.testing.assert_array_almost_equal(siso_sim.B, siso_ini.B)
236+
np.testing.assert_array_almost_equal(siso_sim.C, siso_ini.C)
237+
np.testing.assert_array_almost_equal(siso_sim.D, siso_ini.D)
238+
239+
# Multi-input, multi-output systems
240+
mimo_ini = ss(
241+
[[-1, 1, 0, 0], [0, -2, 1, 0], [0, 0, -3, 1], [0, 0, 0, -4]],
242+
[[1, 0], [0, 0], [0, 1], [1, 1]],
243+
[[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]],
244+
np.zeros((3, 2)))
245+
246+
# Simple transformation: row/col flips + scaling
247+
mimo_txf = np.array(
248+
[[0, 1, 0, 0], [2, 0, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]])
249+
250+
# Transform the system and transform it back
251+
mimo_sim = similarity_transform(mimo_ini, mimo_txf)
252+
mimo_new = similarity_transform(mimo_sim, np.linalg.inv(mimo_txf))
253+
np.testing.assert_array_almost_equal(mimo_new.A, mimo_ini.A)
254+
np.testing.assert_array_almost_equal(mimo_new.B, mimo_ini.B)
255+
np.testing.assert_array_almost_equal(mimo_new.C, mimo_ini.C)
256+
np.testing.assert_array_almost_equal(mimo_new.D, mimo_ini.D)
257+
258+
# Make sure rescaling by identify does nothing
259+
mimo_new = similarity_transform(mimo_ini, np.eye(4))
260+
np.testing.assert_array_almost_equal(mimo_new.A, mimo_ini.A)
261+
np.testing.assert_array_almost_equal(mimo_new.B, mimo_ini.B)
262+
np.testing.assert_array_almost_equal(mimo_new.C, mimo_ini.C)
263+
np.testing.assert_array_almost_equal(mimo_new.D, mimo_ini.D)
264+
265+
# Time rescaling
266+
mimo_tim = similarity_transform(mimo_ini, np.eye(4), timescale=0.3)
267+
mimo_new = similarity_transform(mimo_tim, np.eye(4), timescale=1/0.3)
268+
np.testing.assert_array_almost_equal(mimo_new.A, mimo_ini.A)
269+
np.testing.assert_array_almost_equal(mimo_new.B, mimo_ini.B)
270+
np.testing.assert_array_almost_equal(mimo_new.C, mimo_ini.C)
271+
np.testing.assert_array_almost_equal(mimo_new.D, mimo_ini.D)
272+
273+
# Time + transformation, in one step
274+
mimo_sim = similarity_transform(mimo_ini, mimo_txf, timescale=0.3)
275+
mimo_new = similarity_transform(mimo_sim, np.linalg.inv(mimo_txf),
276+
timescale=1/0.3)
277+
np.testing.assert_array_almost_equal(mimo_new.A, mimo_ini.A)
278+
np.testing.assert_array_almost_equal(mimo_new.B, mimo_ini.B)
279+
np.testing.assert_array_almost_equal(mimo_new.C, mimo_ini.C)
280+
np.testing.assert_array_almost_equal(mimo_new.D, mimo_ini.D)
281+
282+
# Time + transformation, in two steps
283+
mimo_sim = similarity_transform(mimo_ini, mimo_txf, timescale=0.3)
284+
mimo_tim = similarity_transform(mimo_sim, np.eye(4), timescale=1/0.3)
285+
mimo_new = similarity_transform(mimo_tim, np.linalg.inv(mimo_txf))
286+
np.testing.assert_array_almost_equal(mimo_new.A, mimo_ini.A)
287+
np.testing.assert_array_almost_equal(mimo_new.B, mimo_ini.B)
288+
np.testing.assert_array_almost_equal(mimo_new.C, mimo_ini.C)
289+
np.testing.assert_array_almost_equal(mimo_new.D, mimo_ini.D)
290+
221291
def suite():
222292
return unittest.TestLoader().loadTestsFromTestCase(TestFeedback)
223293

0 commit comments

Comments
 (0)
0