|
7 | 7 | importpytest |
8 | 8 |
|
9 | 9 |
|
10 | | -fromcontrolimportStateSpace,forced_response,tf,rss,c2d |
| 10 | +fromcontrolimportStateSpace,impulse_response,forced_response,tf,rss,c2d |
11 | 11 | fromcontrol.exceptionimportControlMIMONotImplemented |
12 | 12 | fromcontrol.tests.conftestimportslycotonly |
13 | | -fromcontrol.modelsimpimportbalred,hsvd,markov,modred |
| 13 | +fromcontrol.modelsimpimportbalred,hsvd,markov,okid,modred |
14 | 14 |
|
15 | 15 |
|
16 | 16 | classTestModelsimp: |
@@ -111,6 +111,142 @@ def testMarkovResults(self, k, m, n): |
111 | 111 | # for k=5, m=n=10: 0.015 % |
112 | 112 | np.testing.assert_allclose(Mtrue,Mcomp,rtol=1e-6,atol=1e-8) |
113 | 113 |
|
| 114 | +deftestOKIDSignature(self): |
| 115 | + |
| 116 | +# Example 6.1, Applied System Identification |
| 117 | +m1,k1,c1=1.,1.,0.01 |
| 118 | +A=np.array([ |
| 119 | + [0.,1.], |
| 120 | + [-k1/m1,-c1/m1], |
| 121 | + ]) |
| 122 | +B=np.array([[0.],[1./m1]]) |
| 123 | +C=np.array([[-k1/m1,-c1/m1]]) |
| 124 | +D=np.array([[1.]]) |
| 125 | +sys=StateSpace(A,B,C,D) |
| 126 | +dt=0.1 |
| 127 | +sysd=sys.sample(dt,method='zoh') |
| 128 | + |
| 129 | +T=np.arange(0,200,dt) |
| 130 | +U=np.random.randn(sysd.B.shape[-1],len(T)) |
| 131 | +response=forced_response(sysd,U=U) |
| 132 | +Y=response.outputs |
| 133 | + |
| 134 | +m=5 |
| 135 | +ir_true=impulse_response(sysd,T=T) |
| 136 | +Htrue=ir_true.outputs[:m+1]*dt |
| 137 | +H=okid(Y,U,m,dt=True) |
| 138 | + |
| 139 | +np.testing.assert_allclose(Htrue,H,atol=1e-1) |
| 140 | + |
| 141 | +# Test mimo example |
| 142 | +# Mechanical Vibrations: Theory and Application, SI Edition, 1st ed. |
| 143 | +# Figure 6.5 / Example 6.7 |
| 144 | +m1,k1,c1=1.,4.,1. |
| 145 | +m2,k2,c2=2.,2.,1. |
| 146 | +k3,c3=6.,2. |
| 147 | + |
| 148 | +A=np.array([ |
| 149 | + [0.,0.,1.,0.], |
| 150 | + [0.,0.,0.,1.], |
| 151 | + [-(k1+k2)/m1, (k2)/m1,-(c1+c2)/m1,c2/m1], |
| 152 | + [(k2)/m2,-(k2+k3)/m2,c2/m2,-(c2+c3)/m2] |
| 153 | + ]) |
| 154 | +B=np.array([[0.,0.],[0.,0.],[1/m1,0.],[0.,1/m2]]) |
| 155 | +C=np.array([[1.0,0.0,0.0,0.0],[0.0,1.0,0.0,0.0]]) |
| 156 | +D=np.zeros((2,2)) |
| 157 | + |
| 158 | +sys=StateSpace(A,B,C,D) |
| 159 | +dt=0.25 |
| 160 | +sysd=sys.sample(dt,method='zoh') |
| 161 | + |
| 162 | +T=np.arange(0,100,dt) |
| 163 | +U=np.random.randn(sysd.B.shape[-1],len(T)) |
| 164 | +response=forced_response(sysd,U=U) |
| 165 | +Y=response.outputs |
| 166 | + |
| 167 | +m=100 |
| 168 | +_,Htrue=impulse_response(sysd,T=dt*(m)) |
| 169 | + |
| 170 | + |
| 171 | +# test array_like |
| 172 | +atol=1e-1 |
| 173 | +H=okid(Y,U,m,dt=dt) |
| 174 | +np.testing.assert_allclose(H,Htrue,atol=atol) |
| 175 | + |
| 176 | +# test array_like, truncate |
| 177 | +H=okid(Y,U,m,dt=dt,truncate=True) |
| 178 | +np.testing.assert_allclose(H,Htrue,atol=atol) |
| 179 | + |
| 180 | +# test array_like, transpose |
| 181 | +HT=okid(Y.T,U.T,m,dt=dt,transpose=True) |
| 182 | +np.testing.assert_allclose(HT,np.transpose(Htrue),atol=atol) |
| 183 | + |
| 184 | +# test response data |
| 185 | +H=okid(response,m,dt=dt) |
| 186 | +np.testing.assert_allclose(H,Htrue,atol=atol) |
| 187 | + |
| 188 | +# test response data |
| 189 | +H=okid(response,m,dt=dt,truncate=True) |
| 190 | +np.testing.assert_allclose(H,Htrue,atol=atol) |
| 191 | + |
| 192 | +# test response data, transpose |
| 193 | +response.transpose=True |
| 194 | +HT=okid(response,m,dt=dt) |
| 195 | +np.testing.assert_allclose(HT,np.transpose(Htrue),atol=atol) |
| 196 | + |
| 197 | +# Make sure markov() returns the right answer |
| 198 | +@pytest.mark.parametrize("k, m, n", |
| 199 | + [(2,2,2), |
| 200 | + (2,5,5), |
| 201 | + (5,2,2), |
| 202 | + (5,5,5), |
| 203 | + (5,10,10)]) |
| 204 | +deftestOkidResults(self,k,m,n): |
| 205 | +# |
| 206 | +# Test over a range of parameters |
| 207 | +# |
| 208 | +# k = order of the system |
| 209 | +# m = number of Markov parameters |
| 210 | +# n = size of the data vector |
| 211 | +# |
| 212 | +# Values *should* match exactly for n = m, otherewise you get a |
| 213 | +# close match but errors due to the assumption that C A^k B = |
| 214 | +# 0 for k > m-2 (see modelsimp.py). |
| 215 | +# |
| 216 | + |
| 217 | +# Generate stable continuous time system |
| 218 | +Hc=rss(k,1,1) |
| 219 | + |
| 220 | +# Choose sampling time based on fastest time constant / 10 |
| 221 | +w,_=np.linalg.eig(Hc.A) |
| 222 | +Ts=np.min(-np.real(w))/10. |
| 223 | + |
| 224 | +# Convert to a discrete time system via sampling |
| 225 | +Hd=c2d(Hc,Ts,'zoh') |
| 226 | + |
| 227 | +# Compute the Markov parameters from state space |
| 228 | +Mtrue=np.hstack([Hd.D]+ [ |
| 229 | +Hd.C @np.linalg.matrix_power(Hd.A,i) @Hd.B |
| 230 | +foriinrange(m-1)]) |
| 231 | + |
| 232 | +Mtrue=np.squeeze(Mtrue) |
| 233 | + |
| 234 | +# Generate input/output data |
| 235 | +T=np.array(range(n))*Ts |
| 236 | +U=np.cos(T)+np.sin(T/np.pi) |
| 237 | +_,Y=forced_response(Hd,T,U,squeeze=True) |
| 238 | +Mcomp=okid(Y,U,m-1) |
| 239 | + |
| 240 | + |
| 241 | +# TODO: Check tol |
| 242 | +# Compare to results from markov() |
| 243 | +# experimentally determined probability to get non matching results |
| 244 | +# with rtot=1e-6 and atol=1e-8 due to numerical errors |
| 245 | +# for k=5, m=n=10: 0.015 % |
| 246 | +np.testing.assert_allclose(Mtrue,Mcomp,atol=1e-1) |
| 247 | + |
| 248 | + |
| 249 | + |
114 | 250 | deftestModredMatchDC(self): |
115 | 251 | #balanced realization computed in matlab for the transfer function: |
116 | 252 | # num = [1 11 45 32], den = [1 15 60 200 60] |
|