Package controlsystems :: Module discretization
[hide private]

Source Code for Module controlsystems.discretization

  1  """Transfer Functions discretization 
  2   
  3  This module implements some numerical methods to discretize the Transfer 
  4  Functions on the time domain. 
  5   
  6  """ 
  7   
  8  __all__ = ['Euler', 'RungeKutta2', 'RungeKutta3', 'RungeKutta4'] 
  9   
 10  #TODO: discretize State-Space models too. 
 11  #TODO: implement more numerical methods 
 12   
 13  from types import Matrix, ZerosMatrix, IdentityMatrix, \ 
 14                    TransferFunction, StateSpace 
 15  from error import ControlSystemsError 
 16   
17 -def Euler(g, sample_time, total_time):
18 """Euler Method 19 20 Returns the points of the step response of the transfer function 'g', 21 discretized with the Euler method, using the sample time 'sample_time' 22 on 'total_time' seconds. For example: 23 24 >>> g = TransferFunction([1], [1, 2, 3]) 25 >>> t, y = Euler(g, 0.01, 10) 26 >>> print t 27 (prints a vector of times 0-10s, with the sample time 0.01s) 28 >>> print y 29 (prints a vector of points, with the same size of 't') 30 31 """ 32 33 if not isinstance(g, TransferFunction): 34 raise ControlSystemsError('Parameter must be a Transfer Fcn.') 35 36 ss = StateSpace(g) 37 38 samples = int(total_time/sample_time) + 1 39 40 t = [sample_time * a for a in range(samples)] 41 42 x = ZerosMatrix(ss.a.rows, 1) 43 y = [] 44 45 eye = IdentityMatrix(ss.a.rows) 46 47 a1 = eye + ss.a.mult(sample_time) 48 a2 = ss.b.mult(sample_time) 49 50 for i in range(samples): 51 x = a1*x + a2 52 y.append((ss.c*x)[0][0] + ss.d[0][0]) 53 54 return t, y
55 56
57 -def RungeKutta2(g, sample_time, total_time):
58 """RungeKutta2 Method 59 60 Returns the points of the step response to the transfer function 'g', 61 discretized with the Runge Kutta (order 2) method, using the sample 62 time 'sample_time' on 'total_time' seconds. For example: 63 64 >>> g = TransferFunction([1], [1, 2, 3]) 65 >>> t, y = RungeKutta2(g, 0.01, 10) 66 >>> print t 67 (prints a vector of times 0-10s, with the sample time 0.01s) 68 >>> print y 69 (prints a vector of points, with the same size of 't') 70 71 """ 72 73 if not isinstance(g, TransferFunction): 74 raise ControlSystemsError('Parameter must be a Transfer Fcn.') 75 76 ss = StateSpace(g) 77 78 samples = int(total_time/sample_time) + 1 79 80 t = [sample_time * a for a in range(samples)] 81 82 x = ZerosMatrix(ss.a.rows, 1) 83 y = [] 84 85 eye = IdentityMatrix(ss.a.rows) 86 87 a1 = ss.a * ss.a 88 a2 = ss.a.mult(2) + a1.mult(sample_time) 89 a3 = a2.mult(0.5) 90 a4 = ss.a.mult(sample_time) 91 a5 = a4*ss.b + ss.b.mult(2) 92 a6 = eye + a3.mult(sample_time) 93 a7 = a5.mult(sample_time/2) 94 95 for i in range(samples): 96 x = a6*x + a7 97 y.append((ss.c*x)[0][0] + ss.d[0][0]) 98 99 return t, y
100 101
102 -def RungeKutta3(g, sample_time, total_time):
103 """RungeKutta3 Method 104 105 Returns the points of the step response to the transfer function 'g', 106 discretized with the Runge Kutta (order 3) method, using the sample 107 time 'sample_time' on 'total_time' seconds. For example: 108 109 >>> g = TransferFunction([1], [1, 2, 3]) 110 >>> t, y = RungeKutta3(g, 0.01, 10) 111 >>> print t 112 (prints a vector of times 0-10s, with the sample time 0.01s) 113 >>> print y 114 (prints a vector of points, with the same size of 't') 115 116 """ 117 118 if not isinstance(g, TransferFunction): 119 raise ControlSystemsError('Parameter must be a Transfer Fcn.') 120 121 ss = StateSpace(g) 122 123 samples = int(total_time/sample_time) + 1 124 125 t = [sample_time * a for a in range(samples)] 126 127 x = ZerosMatrix(ss.a.rows, 1) 128 y = [] 129 130 eye = IdentityMatrix(ss.a.rows) 131 132 a1 = ss.a.mult(sample_time) # A*T 133 a2 = ss.b.mult(sample_time) # B*T 134 a3 = (ss.a * ss.a).mult(sample_time * sample_time) # A^2*T^2 135 a4 = a3.mult(0.5) # (A^2*T^2)/2 136 a5 = (ss.a * ss.b).mult(sample_time * sample_time) # A*B*T^2 137 a6 = a5.mult(0.5) # (A*B*T^2)/2 138 a7 = a3.mult(3.0/4.0) # (A^2*T^2)*(3/4) 139 a8 = (ss.a*ss.a*ss.a).mult(sample_time * sample_time * sample_time) 140 a9 = a8.mult(3.0/8.0) 141 a10 = (ss.a*ss.a*ss.b).mult(sample_time * sample_time * sample_time) 142 a11 = a10.mult(3.0/8.0) 143 a12 = a5.mult(3.0/4.0) + a2 144 145 for i in range(samples): 146 k1 = a1*x + a2 147 k2 = (a1 + a4)*x + a6 + a2 148 k3 = (a1 + a7 + a9)*x + a11 + a12 149 150 x = x + (k1.mult(2.0/9.0) + k2.mult(1.0/3.0) + k3.mult(4.0/9.0)) 151 152 y.append((ss.c*x)[0][0] + ss.d[0][0]) 153 154 return t, y
155 156
157 -def RungeKutta4(g, sample_time, total_time):
158 """RungeKutta4 Method 159 160 Returns the points of the step response to the transfer function 'g', 161 discretized with the Runge Kutta (order 4) method, using the sample 162 time 'sample_time' on 'total_time' seconds. For example: 163 164 >>> g = TransferFunction([1], [1, 2, 3]) 165 >>> t, y = RungeKutta4(g, 0.01, 10) 166 >>> print t 167 (prints a vector of times 0-10s, with the sample time 0.01s) 168 >>> print y 169 (prints a vector of points, with the same size of 't') 170 171 """ 172 173 if not isinstance(g, TransferFunction): 174 raise ControlSystemsError('Parameter must be a Transfer Fcn.') 175 176 ss = StateSpace(g) 177 178 samples = int(total_time/sample_time) + 1 179 180 t = [sample_time * a for a in range(samples)] 181 182 x = ZerosMatrix(ss.a.rows, 1) 183 y = [] 184 185 eye = IdentityMatrix(ss.a.rows) 186 187 a1 = ss.a.mult(sample_time) # A*T 188 a2 = ss.b.mult(sample_time) # B*T 189 a3 = (ss.a * ss.a).mult(sample_time * sample_time) # A^2*T^2 190 a4 = a3.mult(0.5) # (A^2*T^2)/2 191 a5 = (ss.a * ss.b).mult(sample_time * sample_time) # A*B*T^2 192 a6 = a5.mult(0.5) # (A*B*T^2)/2 193 a7 = (ss.a*ss.a*ss.a).mult(sample_time * sample_time * sample_time) 194 a8 = a7.mult(0.25) 195 a9 = (ss.a*ss.a*ss.b).mult(sample_time * sample_time * sample_time) 196 a10 = a9.mult(0.25) 197 a11 = a7.mult(0.5) 198 a12 = (ss.a*ss.a*ss.a*ss.a).mult(sample_time * sample_time * \ 199 sample_time * sample_time) 200 a13 = a12.mult(0.25) 201 a14 = ss.a.mult(sample_time * sample_time) 202 a15 = a14*ss.b 203 a16 = a9.mult(0.5) 204 a17 = (ss.a*ss.a*ss.a*ss.b).mult(sample_time * sample_time * \ 205 sample_time * sample_time) 206 a18 = a17.mult(0.25) 207 208 for i in range(samples): 209 k1 = a1*x + a2 210 k2 = (a1 + a4)*x + a6 + a2 211 k3 = (a1 + a4 + a8)*x + a2 + a6 + a10; 212 k4 = (a1 + a3 + a11 + a13)*x + a15 +a16 + a18 + a2 213 214 x = x + k1.mult(1.0/6.0) + k2.mult(1.0/3.0) + k3.mult(1.0/3.0) \ 215 + k4.mult(1.0/6.0) 216 217 y.append((ss.c*x)[0][0] + ss.d[0][0]) 218 219 return t, y
220