test_sidles_sigg¶

wield.control.test.test_sidles_sigg

Examples demonstrating the use of wield.control objects to calculate Sidles-Sigg instabilities in a number of ways.

First, starting with the state space model of the free torsional spring and treating the hard and soft mode separately 1) Modifying the free state space model 2) Closing the radiation pressure feedback loop around the free plant 3) Using AAA to recover a zpk or state space representation from the frequency response of an absurdly undersampled result of 1)

Second, treating the ITM and ETM together as a MIMO system at the same time and 4) Modifying the MIMO statespace with the undiagonalized stiffness matrix 5) Applying the radiation pressure feedback through a MIMO state space feedback D matrix

Control feedback can be applied in the same way as the radiation pressure feedback in this example is.

Provided by Kevin Kuns

pytest-html report

Functions

test_hard_soft_modes(capture)

This is a pytest needing documentation

Details

test_hard_soft_modes(capture)[source][github]¶

This is a pytest needing documentation

code
  1def test_hard_soft_modes(capture):
  2    F_Hz = np.logspace(-2, 1, 1000)
  3
  4    # Parameters
  5    F0_Hz = 0.5
  6    w0_rad_s = 2 * np.pi * F0_Hz
  7    Q = 100
  8    I_kgm2 = 2.73
  9    Parm_W = 2e6
 10    Larm_m = 4e3
 11    Ri_m = 1934
 12    Re_m = 2245
 13    gi = 1 - Larm_m / Ri_m
 14    ge = 1 - Larm_m / Re_m
 15
 16    # Hard and soft torsional stiffness
 17    k0 = 2 * Parm_W * Larm_m / (scc.c * (gi * ge - 1))
 18    kh = k0 * (ge + gi - np.sqrt((ge - gi)**2 + 4)) / 2
 19    ks = k0 * (ge + gi + np.sqrt((ge - gi)**2 + 4)) / 2
 20
 21    ###########################################################################
 22    # SISO state space
 23    ###########################################################################
 24
 25    # Free suspension state space matricies
 26    Afree_siso = np.array([
 27        [0, 1],
 28        [-w0_rad_s**2, -w0_rad_s / Q],
 29    ])
 30    B_siso = np.array([
 31        [0],
 32        [1 / I_kgm2],
 33    ])
 34    C_siso = np.array([[1, 0]])
 35    D_siso = np.array([[0]])
 36
 37    # Hard and soft modes shift the free torsional stiffness
 38    Khard = np.array([
 39        [0, 0],
 40        [kh / I_kgm2, 0],
 41    ])
 42    Ksoft = np.array([
 43        [0, 0],
 44        [ks / I_kgm2, 0],
 45    ])
 46    Ahard = Afree_siso - Khard
 47    Asoft = Afree_siso - Ksoft
 48
 49    # Make SISOStateSpace objects using the modified state space matricies
 50    free_siso = SISO.SISOStateSpace(Afree_siso, B_siso, C_siso, D_siso)
 51    hard_siso = SISO.SISOStateSpace(Ahard, B_siso, C_siso, D_siso)
 52    soft_siso = SISO.SISOStateSpace(Asoft, B_siso, C_siso, D_siso)
 53
 54    # Find the hard and soft mode poles
 55    # To find the poles, SISOStateSpace objects can be converted into ZPK
 56    # objects using siso.asZPK. Likewise, ZPK objects can be converted into
 57    # SISOStateSpace objects using zpk.asSS
 58    dprint('Hard poles [Hz]', hard_siso.asZPK.p / (2 * np.pi))
 59    dprint('Soft poles [Hz]', soft_siso.asZPK.p / (2 * np.pi))
 60    dprint('Soft stable?', np.all(soft_siso.asZPK.p.real < 0))
 61
 62    # SISO zpk objects can also be defined directly
 63    # frequencies are given in rad/s by default but can be given in Hz with
 64    # the angular=False keyword
 65    Fp_Hz = np.array([
 66        -F0_Hz / (2 * Q) * (1 + 1j * np.sqrt(4 * Q**2 - 1)),
 67        -F0_Hz / (2 * Q) * (1 - 1j * np.sqrt(4 * Q**2 - 1))
 68    ])
 69    free_zpk = SISO.zpk([], Fp_Hz, 1 / I_kgm2, angular=False)
 70    np.testing.assert_almost_equal(
 71        np.sort(free_zpk.p), np.sort(free_siso.asZPK.p))
 72
 73    # The hard and soft mode plants can also be found by closing the radiation
 74    # pressure loop: the OLG is -k * (free plant)
 75    # SISO objects can be added, multiplied, and divided like variables
 76    hard_rp_loop = (1 / (1 + kh * free_siso)) * free_siso
 77    soft_rp_loop = (1 / (1 + ks * free_siso)) * free_siso
 78
 79    # Plot comparisons
 80    # Frequency response is calculated with a fresponse object. The frequency
 81    # vector can be specified either in Hz (with the f keyword), in rad/s (with
 82    # the w keyword), or in the s-domain (with the s keyword).
 83    # The complex numerical array is given by the tf attribute
 84    fig = plotTF(
 85        F_Hz, free_siso.fresponse(f=F_Hz).tf,
 86        label='Free', c='xkcd:kelly green',
 87    )
 88    plotTF(
 89        F_Hz, hard_siso.fresponse(f=F_Hz).tf, *fig.axes,
 90        label='Hard (SISO statespace)', c='xkcd:tangerine',
 91    )
 92    plotTF(
 93        F_Hz, soft_siso.fresponse(f=F_Hz).tf, *fig.axes,
 94        label='Soft (SISO statespace)', c='xkcd:burgundy',
 95    )
 96    plotTF(
 97        F_Hz, hard_rp_loop.fresponse(f=F_Hz).tf, *fig.axes,
 98        label='Hard (close RP loop)',
 99        ls='--', c='xkcd:royal purple',
100    )
101    plotTF(
102        F_Hz, soft_rp_loop.fresponse(f=F_Hz).tf, *fig.axes,
103        label='Soft (close RP loop)',
104        ls='--', c='xkcd:sky blue',
105    )
106    fig.axes[1].legend(loc='upper left')
107    fig.set_size_inches((6, 6.4))
108    fig.savefig(tjoin('compare_siso.pdf'))
109
110    ###########################################################################
111    # AAA fit
112    ###########################################################################
113
114    # Generate some absurdly undersampled data to fit
115    F_fit_Hz = np.logspace(-1, 1, 5)
116    hard_fit_data = hard_siso.fresponse(f=F_fit_Hz).tf
117    soft_fit_data = soft_siso.fresponse(f=F_fit_Hz).tf
118
119    # Use AAA to fit this data
120    hard_fit = tfAAA(F_fit_Hz, hard_fit_data)
121    soft_fit = tfAAA(F_fit_Hz, soft_fit_data)
122
123    # Define a zpk object from these fits
124    # AAA uses the IIRrational normalization, which needs to be specified
125    # The standard normalization ('scipy') is the default
126    hard_aaa = SISO.zpk(*hard_fit.zpk, convention='iirrational')
127    soft_aaa = SISO.zpk(*soft_fit.zpk, convention='iirrational')
128
129    # Check that the fit recovered the poles of the orginal SISO model
130    np.testing.assert_almost_equal(
131        np.sort(hard_siso.asZPK.p), np.sort(hard_aaa.p))
132    np.testing.assert_almost_equal(
133        np.sort(soft_siso.asZPK.p), np.sort(soft_aaa.p))
134
135    # Plot the data on top of the fits
136    fig = plotTF(
137        F_Hz, hard_aaa.fresponse(f=F_Hz).tf,
138        label='Hard AAA fit', c='xkcd:cerulean',
139    )
140    plotTF(
141        F_Hz, soft_aaa.fresponse(f=F_Hz).tf, *fig.axes,
142        label='Soft AAA fit', c='xkcd:tangerine',
143    )
144    plotTF(
145        F_fit_Hz, hard_fit_data, *fig.axes,
146        label='Hard AAA fit data',
147        ls='', marker='o', c='xkcd:cerulean')
148    plotTF(
149        F_fit_Hz, soft_fit_data, *fig.axes,
150        label='Soft AAA fit data',
151        ls='', marker='o', c='xkcd:tangerine',
152    )
153    fig.axes[1].legend(loc='upper left')
154    fig.axes[0].set_xlim(F_Hz[0], F_Hz[-1])
155    fig.set_size_inches((6, 6.4))
156    fig.savefig(tjoin('compare_aaa.pdf'))
157
158    ###########################################################################
159    # MIMO state space
160    ###########################################################################
161
162    # Free suspension state space matrices
163    # State variables are
164    # (ETM angle, ITM angle, ETM angular velocity, ITM angular velocity)
165    eye = np.eye(2)
166    Afree_mimo = np.block([
167        [0 * eye, eye],
168        [-w0_rad_s**2 * eye, -w0_rad_s / Q * eye],
169    ])
170    B_mimo = np.block([
171        [0 * eye],
172        [1 / I_kgm2 * eye],
173    ])
174    C_mimo = np.block([[eye, 0 * eye]])
175    D_mimo = 0 * eye
176
177    # RP torsional stiffness matrix
178    Krp = k0 * np.array([
179        [gi, 1],
180        [1, ge],
181    ])
182    K_mimo = np.block([
183        [0 * eye, 0 * eye],
184        [Krp / I_kgm2, 0 * eye],
185    ])
186    A_mimo = Afree_mimo - K_mimo
187
188    # Make MIMOStateSpace objects
189    # Input/output degrees of freedom/test points are defined by either lists
190    # or dictionaries specifying the row/column indices to which they
191    # correspond. They can have different names but are the same in this example.
192    dofs = {'etm': 0, 'itm': 1}
193    # dofs = ['etm', 'itm']
194    free_mimo = MIMO.MIMOStateSpace(
195        Afree_mimo, B_mimo, C_mimo, D_mimo,
196        inputs=dofs,
197        outputs=dofs,
198    )
199    ss_mimo = MIMO.MIMOStateSpace(
200        A_mimo, B_mimo, C_mimo, D_mimo,
201        inputs=dofs,
202        outputs=dofs,
203    )
204
205    # The radiation pressure modified state space can also be defined by
206    # specifying feedback connections in the free model to define a feedback
207    # D matrix
208    connections = {
209        ('etm', 'etm'): -Krp[0, 0],
210        ('etm', 'itm'): -Krp[0, 1],
211        ('itm', 'etm'): -Krp[1, 0],
212        ('itm', 'itm'): -Krp[1, 1],
213    }
214    fback_ss = free_mimo.feedback_connect(connections=connections)
215
216    # These MIMOStateSpace objects are in the mirror basis. (Choosing different
217    # B and C matrices could have put them in the hard/soft basis instead, in
218    # which case the dofs dictionary should have been changed to have keys
219    # 'hard' and 'soft' instead).
220    # Individual SISOStateSpace objects can be extracted like
221    free_etm_siso = free_mimo.siso('etm', 'etm')
222
223    # # MIMO fresponse (in the mirror basis) is calculated like SISO is.
224    # # SISO fresponse can also be extracted from these MIMO fresponse objects
225    mimo_response = ss_mimo.fresponse(f=F_Hz)
226    fback_response = fback_ss.fresponse(f=F_Hz)
227
228    # Plot results in the mirror basis
229    fig = plotTF(
230        F_Hz, free_etm_siso.fresponse(f=F_Hz).tf,
231        label='Free', c='xkcd:kelly green',
232    )
233    plotTF(
234        F_Hz, mimo_response.siso('etm', 'etm').tf, *fig.axes,
235        label='ETM to ETM', c='xkcd:cerulean',
236    )
237    plotTF(
238        F_Hz, mimo_response.siso('etm', 'itm').tf, *fig.axes,
239        label='ITM to ETM', c='xkcd:tangerine',
240    )
241    fig.axes[1].legend(loc='upper left')
242    fig.set_size_inches((6, 6.4))
243    fig.savefig(tjoin('mimo_mirror_basis.pdf'))
244
245    # Analyze results in the hard/soft basis.
246    # The full numerical MIMO plant (now a matrix) is given by the tf attribute
247    # The hard and soft eigenvectors are
248    _, eigv = np.linalg.eig(Krp)
249    vhard = eigv[:, 0]
250    vsoft = eigv[:, 1]
251
252    # Convert to hard/soft basis
253    def to_hard_soft(fresponse):
254        return Bunch(
255            hard = vhard.T @ fresponse.tf @ vhard,
256            soft = vsoft.T @ fresponse.tf @ vsoft,
257        )
258
259    mimo_plants = to_hard_soft(mimo_response)
260    fback_plants = to_hard_soft(fback_response)
261
262    # Plot comparisons
263    fig = plotTF(
264        F_Hz, free_etm_siso.fresponse(f=F_Hz).tf,
265        label='Free', c='xkcd:kelly green',
266    )
267    plotTF(
268        F_Hz, mimo_plants.hard, *fig.axes,
269        label='Hard (MIMO statespace)', c='xkcd:tangerine',
270    )
271    plotTF(
272        F_Hz, mimo_plants.soft, *fig.axes,
273        label='Soft (MIMO statespace)', c='xkcd:burgundy',
274    )
275    plotTF(
276        F_Hz, fback_plants.hard, *fig.axes,
277        label='Hard (MIMO feedback)',
278        ls='--', c='xkcd:royal purple',
279    )
280    plotTF(
281        F_Hz, fback_plants.soft, *fig.axes,
282        label='Soft (MIMO feedback)',
283        ls='--', c='xkcd:sky blue',
284    )
285    fig.axes[1].legend(loc='upper left')
286    fig.set_size_inches((6, 6.4))
287    fig.savefig(tjoin('compare_mimo.pdf'))
288
289    ###########################################################################
290    # Check that all of these methods are equavalent
291    ###########################################################################
292
293    def assert_equal(hard, soft):
294        np.testing.assert_almost_equal(
295            hard_siso.fresponse(f=F_Hz).tf,
296            hard,
297        )
298        np.testing.assert_almost_equal(
299            soft_siso.fresponse(f=F_Hz).tf,
300            soft,
301        )
302
303    assert_equal(
304        hard_rp_loop.fresponse(f=F_Hz).tf,
305        soft_rp_loop.fresponse(f=F_Hz).tf,
306    )
307    assert_equal(
308        hard_aaa.fresponse(f=F_Hz).tf,
309        soft_aaa.fresponse(f=F_Hz).tf,
310    )
311    assert_equal(
312        mimo_plants.hard,
313        mimo_plants.soft,
314    )
315    assert_equal(
316        fback_plants.hard,
317        fback_plants.soft,
318    )
pytest information

This code is wrapped in a pytest function using conventions detailed in Pytest Conventions. The full name of this test, as known by the documentation, is:

wield.control.test.test_sidles_sigg.test_hard_soft_modes

The full name is useful when building documentation, to link a reference to this page using :func:`name`, or directly include it with an autofunction directive. The collapse nodes below show every instance of the test run. There may only be one, but if the test was run multiple times through pytest parametrizations, the list can be longer.

test_hard_soft_modes
output
status: passed
duration: 4.668s
Captured stderr call
'Hard poles [Hz]' array([-0.0025+2.46194271j, -0.0025-2.46194271j])
'Soft poles [Hz]' array([ 0.07427487, -0.07927487])
'Soft stable?' False