Skip to content

Commit

Permalink
Merge 47e780c into 9d2bfe3
Browse files Browse the repository at this point in the history
  • Loading branch information
@bnavigator
bnavigator committed Apr 10, 2020
2 parents 9d2bfe3 + 47e780c commit 3c197ff
Show file tree
Hide file tree
Showing 4 changed files with 144 additions and 76 deletions.
46 changes: 29 additions & 17 deletions slycot/synthesis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -489,16 +489,22 @@ def sb02mt(n,m,B,R,A=None,Q=None,L=None,fact='N',jobl='Z',uplo='U',ldwork=None):
if fact == 'N':
out = _wrapper.sb02mt_n(n,m,B,R,uplo=uplo,ldwork=ldwork)
if out is None:
raise ValueError('fact must be either C or N.')
e = ValueError('fact must be either C or N.')
e.info = -3
raise e
else:
if A is None or Q is None or L is None:
raise ValueError('matrices A,Q and L are required if jobl is not Z.')
e = ValueError('matrices A,Q and L are required if jobl is not Z.')
e.info = -7
raise e
if fact == 'C':
out = _wrapper.sb02mt_cl(n,m,A,B,Q,R,L,uplo=uplo)
if fact == 'N':
out = _wrapper.sb02mt_nl(n,m,A,B,Q,R,L,uplo=uplo,ldwork=ldwork)
if out is None:
raise ValueError('fact must be either C or N.')
e = ValueError('fact must be either C or N.')
e.info = -3
raise e
if out[-1] < 0:
error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1]
e = ValueError(error_text)
Expand Down Expand Up @@ -850,7 +856,9 @@ def sb03md(n,C,A,U,dico,job='X',fact='N',trana='N',ldwork=None):
if ldwork is None:
ldwork = max(2*n*n,3*n)
if dico != 'C' and dico != 'D':
raise ValueError('dico must be either D or C')
e = ValueError('dico must be either D or C')
e.info = -1
raise e
out = _wrapper.sb03md(dico,n,C,A,U,job=job,fact=fact,trana=trana,ldwork=ldwork)
if out[-1] < 0:
error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1]
Expand Down Expand Up @@ -1041,7 +1049,9 @@ def sb03od(n,m,A,Q,B,dico,fact='N',trans='N',ldwork=None):
elif m == 0:
ldwork = 1
if dico != 'C' and dico != 'D':
raise ValueError('dico must be either D or C')
e = ValueError('dico must be either D or C')
e.info = -1
raise e
out = _wrapper.sb03od(dico,n,m,A,Q,B,fact=fact,trans=trans,ldwork=ldwork)
if out[-1] < 0:
error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1]
Expand Down Expand Up @@ -1658,7 +1668,7 @@ def sb10dd(n,m,np,ncon,nmeas,gamma,A,B,C,D,tol=0.0,ldwork=None):
LW4 = 13*n*n + m*m + (8*n+m+m2+2*np2)*(m2+np2) + 6*n + n*(m+np2) + max(14*n+23,16*n,2*n+m2+np2,3*(m2+np2))
ldwork = max(LW1,LW2,LW3,LW4)
out = _wrapper.sb10dd(n,m,np,ncon,nmeas,gamma,A,B,C,D,tol,ldwork)

if out[-1] != 0:
if out[-1] < 0:
error_text = "The following argument had an illegal value: "\
Expand Down Expand Up @@ -1837,25 +1847,25 @@ def sb10hd(n,m,np,ncon,nmeas,A,B,C,D,tol=0.0,ldwork=None):
raise e

return out[:-1]

def sb10jd(n,m,np,A,B,C,D,E,ldwork=None):
""" A,B,C,D = sb10jd(n,m,np,A,B,C,D,E,[ldwork])
To convert the descriptor state-space system
E*dx/dt = A*x + B*u
y = C*x + D*u
into regular state-space form
dx/dt = Ad*x + Bd*u
y = Cd*x + Dd*u .
Required arguments:
n : input int
The order of the descriptor system. n >= 0.
m : input int
The column size of the matrix B. m >= 0.
The column size of the matrix B. m >= 0.
np : input int
The row size of the matrix C. np >= 0.
A : rank-2 array('d') with bounds (n,n)
Expand Down Expand Up @@ -1890,7 +1900,7 @@ def sb10jd(n,m,np,A,B,C,D,E,ldwork=None):
contains the output matrix Cd of the converted system.
D : rank-2 array('d') with bounds (np,m)
The leading NP-by-M part of this array contains
the matrix Dd of the converted system.
the matrix Dd of the converted system.
"""

hidden = ' (hidden by the wrapper)'
Expand All @@ -1900,7 +1910,7 @@ def sb10jd(n,m,np,A,B,C,D,E,ldwork=None):
ldwork = max(1, 2 * n * n + 2 * n + n * max(5, n + m + np))

A,B,C,D,nsys,info = _wrapper.sb10jd(n,m,np,A,B,C,D,E,ldwork)

if info < 0:
error_text = "The following argument had an illegal value: "+arg_list[-info-1]
e = ValueError(error_text)
Expand Down Expand Up @@ -2694,7 +2704,7 @@ def sg03bd(n,m,A,E,Q,Z,B,dico,fact='N',trans='N',ldwork=None):
Return objects
______________
U : rank-2 array('d'), shape (n,n)
The leading n-by-b part of this array contains
the Cholesky factor U of the solution matrix X of the
Expand Down Expand Up @@ -2746,7 +2756,9 @@ def sg03bd(n,m,A,E,Q,Z,B,dico,fact='N',trans='N',ldwork=None):
if ldwork is None:
ldwork = max(1,4*n,6*n-6)
if dico != 'C' and dico != 'D':
raise ValueError('dico must be either D or C')
e = ValueError('dico must be either D or C')
e.info = -1
raise e
out = _wrapper.sg03bd(dico,n,m,A,E,Q,Z,B,fact=fact,trans=trans,ldwork=ldwork)
if out[-1] < 0:
error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1]
Expand Down
14 changes: 0 additions & 14 deletions slycot/tests/test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,20 +43,6 @@ def test_sb02ad(self):
self.assertAlmostEqual(Ac[1][0], 1)
self.assertAlmostEqual(Ac[1][1], -3)

def test_td04ad_static(self):
"""Regression: td04ad (TFM -> SS transformation) for static TFM"""
import numpy as np
from itertools import product
# 'C' fails on static TFs
for nout,nin,rc in product(range(1,6),range(1,6),['R']):
num = np.reshape(np.arange(nout*nin),(nout,nin,1))
if rc == 'R':
den = np.reshape(np.arange(1,1+nout),(nout,1))
else:
den = np.reshape(np.arange(1,1+nin),(nin,1))
index = np.tile([0],den.shape[0])
nr,a,b,c,d = transform.td04ad(rc,nin,nout,index,den,num)


if __name__ == "__main__":
unittest.main()
148 changes: 109 additions & 39 deletions slycot/tests/test_td04ad.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,19 +3,18 @@
# test_td04ad.py - test suite for tf -> ss conversion
# RvP, 04 Jun 2018

from __future__ import print_function
from __future__ import print_function, division

import unittest
from slycot import transform
import numpy as np

from numpy.testing import assert_raises, assert_almost_equal

class TestTf2SS(unittest.TestCase):

def test_td04ad_case1(self):
"""td04ad: Convert with both 'C' and 'R' options"""
def test_td04ad_c(self):
"""td04ad: Convert with 'C' option"""

# for octave:
"""
num = { [0.0, 0.0, 1.0 ], [ 1.0, 0.0 ];
Expand All @@ -25,80 +24,155 @@ def test_td04ad_case1(self):
[1.0, 0.4, 3.0], [ 1.0, 1.0 ];
[1.0, 0.4, 3.0], [ 1.0, 1.0 ]};
"""

# common denominators for the inputs
n = 2

m = 2
p = 3
d = 3
num = np.array([
[ [0.0, 0.0, 1.0 ], [ 1.0, 0.0, 0.0 ] ],
[ [3.0, -1.0, 1.0 ], [ 0.0, 1.0, 0.0 ] ],
[ [0.0, 0.0, 1.0], [ 0.0, 2.0, 0.0 ] ] ])
p, m, d = num.shape
[ [0.0, 0.0, 1.0], [1.0, 0.0, 0.0] ],
[ [3.0, -1.0, 1.0], [0.0, 1.0, 0.0] ],
[ [0.0, 0.0, 1.0], [0.0, 2.0, 0.0] ] ])

numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float)
numc[:p,:m,:] = num

denc = np.array(
[ [1.0, 0.4, 3.0], [ 1.0, 1.0, 0.0 ] ])
indc = np.array(
[ 2, 1 ], dtype=int)
denr = np.array(
[ [1.0, 0.4, 3.0], [ 1.0, 1.0, 0.0 ], [1.0, 0.0, 0.0] ])
indr = np.array(
[ 2, 1, 0 ], dtype=int)

n, A, B, C, D = transform.td04ad('C', 2, 3, indc, denc, numc)
nref = 3
Aref = np.array([ [-1, 0, 0],
[ 0, -0.4, -0.3],
[ 0, 10, 0] ])
Bref = np.array([ [0, -1],
[1, 0],
[0, 0] ])
Cref = np.array([ [1, 0, 0.1],
[-1, -2.2, -0.8],
[-2, 0, 0.1] ])
Dref = np.array([ [0, 1],
[3, 0],
[0, 0] ])

nr, A, B, C, D = transform.td04ad('C', m, p, indc, denc, numc)
#print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
Ac = [ [-1, 0, 0], [ 0, -0.4, -0.3], [ 0, 10, 0]]
Bc = [ [0, -1] ,[ 1 , 0], [ 0, 0]]
Cc = [ [1, 0, 0.1], [-1, -2.2, -0.8], [ -2, 0, 0.1] ]
Dc = [ [0, 1], [ 3, 0], [ 0, 0]]
np.testing.assert_array_almost_equal(A, Ac)
np.testing.assert_array_almost_equal(B, Bc)
np.testing.assert_array_almost_equal(C, Cc)
np.testing.assert_array_almost_equal(D, Dc)
np.testing.assert_equal(nref, nr)
# the returned state space representation is not guaranteed to
# be of one form for all architectures, so we transform back
# to tf and check for equality then
_, _, _, _, _, dcoeff, ucoeff = transform.tb04ad(
nr, m, p, A, B, C, D)
_, _, _, _, _, dcoeffref, ucoeffref = transform.tb04ad(
nref, m, p, Aref, Bref, Cref, Dref)
np.testing.assert_array_almost_equal(dcoeff,dcoeffref)
np.testing.assert_array_almost_equal(ucoeff,ucoeffref)


def test_td04ad_r(self):
"""td04ad: Convert with 'R' option"""

resr = transform.td04ad('R', 2, 3, indr, denr, num)
""" example program from
http://slicot.org/objects/software/shared/doc/TD04AD.html"""

m = 2
p = 2
rowcol = 'R'
index = [3, 3]
dcoeff = np.array([ [1.0, 6.0, 11.0, 6.0], [1.0, 6.0, 11.0, 6.0] ])

ucoeff = np.array([ [[1.0, 6.0, 12.0, 7.0], [0.0, 1.0, 4.0, 3.0]],
[[0.0, 0.0, 1.0, 1.0], [1.0, 8.0, 20.0, 15.0]] ])

nref = 3

Aref = np.array([ [ 0.5000, -0.8028, 0.9387],
[ 4.4047, -2.3380, 2.5076],
[-5.5541, 1.6872, -4.1620] ])
Bref = np.array([ [-0.2000, -1.2500],
[ 0.0000, -0.6097],
[ 0.0000, 2.2217] ])
Cref = np.array([ [0.0000, -0.8679, 0.2119],
[0.0000, 0.0000, 0.9002] ])
Dref = np.array([ [1.0000, 0.0000],
[0.0000, 1.0000] ])

nr, A, B, C, D = transform.td04ad(rowcol, m, p, index, dcoeff, ucoeff)
#print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
np.testing.assert_equal(nref, nr)
# order of states is not guaranteed, so we reorder the reference
rindex = np.flip(np.argsort(np.diag(A)))
Arref = Aref[rindex, :][:, rindex]
Brref = Bref[rindex, :]
Crref = Cref[:, rindex]
Drref = Dref
np.testing.assert_array_almost_equal(A, Arref,decimal=4)
np.testing.assert_array_almost_equal(B, Brref,decimal=4)
np.testing.assert_array_almost_equal(C, Crref,decimal=4)
np.testing.assert_array_almost_equal(D, Drref,decimal=4)


def test_staticgain(self):
"""td04ad: Convert a transferfunction to SS with only static gain"""

# 2 inputs, 3 outputs? columns share a denominator
num = np.array([ [ [1.0], [2.0] ],
[ [0.2], [4.3] ],
[ [1.2], [3.2] ] ])
p, m, d = num.shape
numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float)
numc[:p,:m,:] = num

# denc, columns share a common denominator
denc = np.array([ [ 1.0], [0.5] ])
Dc = (num / denc).reshape((3,2))
idxc = np.zeros((2,), dtype=int)

# denr, rows share a common denominator
denr = np.array([ [1.0], [0.5], [3.0] ])
idxr = np.zeros((3,), dtype=int)
Dr = (num / denr[:, np.newaxis]).reshape((3,2))

# fails with:
# On entry to TB01XD parameter number 5 had an illegal value

n, A, B, C, D = transform.td04ad('C', 2, 3, idxc, denc, numc)
#print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
self.assertEqual(A.shape, (0,0))
self.assertEqual(B.shape, (0,2))
self.assertEqual(C.shape, (3,0))
np.testing.assert_array_almost_equal(D, Dc)

n, A, B, C, D = transform.td04ad('R', 2, 3, idxr, denr, num)
#print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
self.assertEqual(A.shape, (0,0))
self.assertEqual(B.shape, (0,2))
self.assertEqual(C.shape, (3,0))
np.testing.assert_array_almost_equal(D, Dr)


def test_td04ad_static(self):
"""Regression: td04ad (TFM -> SS transformation) for static TFM"""
from itertools import product
for nout, nin, rc in product(range(1, 6), range(1, 6), ['R', 'C']):
Dref = np.zeros((nout, nin))
if rc == 'R':
num = np.reshape(np.arange(nout * nin), (nout, nin, 1))
den = np.reshape(np.arange(1, 1 + nout), (nout, 1))
index = np.repeat(0, nout)
Dref = num[:nout, :nin, 0] / np.broadcast_to(den, (nout, nin))
else:
maxn = max(nout, nin)
num = np.zeros((maxn, maxn, 1))
num[:nout, :nin, 0] = np.reshape(
np.arange(nout * nin), (nout, nin))
den = np.reshape(np.arange(1, 1 + nin), (nin, 1))
index = np.repeat(0, nin)
Dref = num[:nout, :nin, 0] / np.broadcast_to(den.T, (nout, nin))
nr, A, B, C, D = transform.td04ad(rc, nin, nout, index, den, num)
np.testing.assert_equal(nr, 0)
for M in [A, B, C]:
np.testing.assert_equal(M, np.zeros_like(M))
np.testing.assert_almost_equal(D, Dref)


def test_mixfeedthrough(self):
"""Test case popping up from control testing"""
Expand All @@ -113,7 +187,7 @@ def test_mixfeedthrough(self):
idxc = np.array([ 1, 0 ])
n, A, B, C, D = transform.td04ad('C', 2, 2, idxc, denc, numc)
np.testing.assert_array_almost_equal(D, np.array([[0, 0],[-0.1, 0]]))

def test_toandfrom(self):

A = np.array([[-3.0]])
Expand All @@ -131,7 +205,7 @@ def test_toandfrom(self):
np.testing.assert_array_almost_equal(A, At)

def test_tfm2ss_6(self):
"""Python version of Fortran test program from
"""Python version of Fortran test program from
-- Bug in TD04AD when ROWCOL='C' #6
This bug was fixed in PR #27"""
m = 1
Expand All @@ -145,10 +219,6 @@ def test_tfm2ss_6(self):
n, A, B, C, D = transform.td04ad('C', m, p, index, dcoeff, ucoeff)
self.assertEqual(n, 0)
np.testing.assert_array_almost_equal(D, np.array([[64]]))

def suite():
return unittest.TestLoader().loadTestsFromTestCase(TestTF2SS)


if __name__ == "__main__":
unittest.main()
Loading

0 comments on commit 3c197ff

Please sign in to comment.