44from .exception import ControlNotImplemented
55from .lti import issiso
66from .statesp import StateSpace
7- from .statefbk import ctrb
7+ from .statefbk import ctrb , obsv
88
99from numpy import zeros ,shape ,poly
1010from numpy .linalg import inv
1111
12- __all__ = ['canonical_form' ,'reachable_form' ]
12+ __all__ = ['canonical_form' ,'reachable_form' , 'observable_form' ]
1313
1414def canonical_form (xsys ,form = 'reachable' ):
1515"""Convert a system into canonical form
@@ -21,7 +21,7 @@ def canonical_form(xsys, form='reachable'):
2121 form : String
2222 Canonical form for transformation. Chosen from:
2323 * 'reachable' - reachable canonical form
24- * 'observable' - observable canonical form [not implemented]
24+ * 'observable' - observable canonical form
2525 * 'modal' - modal canonical form [not implemented]
2626
2727 Returns
@@ -35,6 +35,8 @@ def canonical_form(xsys, form='reachable'):
3535# Call the appropriate tranformation function
3636if form == 'reachable' :
3737return reachable_form (xsys )
38+ elif form == 'observable' :
39+ return observable_form (xsys )
3840else :
3941raise ControlNotImplemented (
4042"Canonical form '%s' not yet implemented" % form )
@@ -85,3 +87,49 @@ def reachable_form(xsys):
8587zsys .C = xsys .C * inv (Tzx )
8688
8789return zsys ,Tzx
90+
91+
92+ def observable_form (xsys ):
93+ """Convert a system into observable canonical form
94+
95+ Parameters
96+ ----------
97+ xsys : StateSpace object
98+ System to be transformed, with state `x`
99+
100+ Returns
101+ -------
102+ zsys : StateSpace object
103+ System in observable canonical form, with state `z`
104+ T : matrix
105+ Coordinate transformation: z = T * x
106+ """
107+ # Check to make sure we have a SISO system
108+ if not issiso (xsys ):
109+ raise ControlNotImplemented (
110+ "Canonical forms for MIMO systems not yet supported" )
111+
112+ # Create a new system, starting with a copy of the old one
113+ zsys = StateSpace (xsys )
114+
115+ # Generate the system matrices for the desired canonical form
116+ zsys .C = zeros (shape (xsys .C ))
117+ zsys .C [0 ,0 ]= 1
118+ zsys .A = zeros (shape (xsys .A ))
119+ Apoly = poly (xsys .A )# characteristic polynomial
120+ for i in range (0 ,xsys .states ):
121+ zsys .A [i ,0 ]= - Apoly [i + 1 ]/ Apoly [0 ]
122+ if (i + 1 < xsys .states ):
123+ zsys .A [i ,i + 1 ]= 1
124+
125+ # Compute the observability matrices for each set of states
126+ Wrx = obsv (xsys .A ,xsys .C )
127+ Wrz = obsv (zsys .A ,zsys .C )
128+
129+ # Transformation from one form to another
130+ Tzx = inv (Wrz )* Wrx
131+
132+ # Finally, compute the output matrix
133+ zsys .B = Tzx * xsys .B
134+
135+ return zsys ,Tzx
