Parameter Estimation of RIG Roll Experiments

Setup and descriptions

  • Without ACM model
  • Turn off wind tunnel
  • Only 1DoF for RIG roll movement
  • Free in any angle and observe
  • Consider RIG roll angle and its derivative as States (in radians) $$X = \begin{pmatrix} \phi_{rig} \\ \dot{\phi}_{rig} \end{pmatrix}$$
  • Observe RIG roll angle and its derivative as Outputs (in degrees) $$Z = \begin{pmatrix} \phi_{rig} \\ \dot{\phi}_{rig} \end{pmatrix}$$
  • Use output error method based on most-likehood(ML) to estimate $$ \theta = \begin{pmatrix} F_{c} \\ f \end{pmatrix} $$

Startup computation engines


In [2]:
%run matt_startup
%run -i matt_utils

button_qtconsole()

#import other needed modules in all used engines
#with dview.sync_imports():
#   import os


importing numpy on engine(s)
importing scipy on engine(s)
importing matplotlib on engine(s)
importing scipy.integrate on engine(s)
importing scipy.interpolate on engine(s)
importing math on engine(s)
importing time on engine(s)
importing itertools on engine(s)
importing random on engine(s)

Data preparation

Load raw data


In [33]:
filename = 'FIWT_Exp014_20150601105647.dat.npz'

def loadData():
    # Read and parse raw data
    global exp_data
    exp_data = np.load(filename)

    # Select colums
    global T_rig, phi_rig
    T_rig = exp_data['data44'][:,0]
    phi_rig = exp_data['data44'][:,2]

loadData()

In [34]:
text_loadData()

Check time sequence and inputs/outputs

Click 'Check data' button to show the raw data.

Click on curves to select time point and push into queue; click 'T/s' text to pop up last point in the queue; and click 'Output' text to print time sequence table.


In [35]:
def checkInputOutputData():

    #check inputs/outputs
    fig, ax = plt.subplots(2,1,True)
    ax[1].plot(T_rig,phi_rig, 'b', picker=2)
    ax[1].set_ylabel('$\phi \/ / \/ ^o/s$')
    ax[1].set_xlabel('$T \/ / \/ s$', picker=True)
    ax[0].set_title('Output', picker=True)
    

    fig.canvas.mpl_connect('pick_event', onPickTime)
    fig.show()
    display(fig)

button_CheckData()

Input $\delta_T$ and focused time ranges

For each section,

  • Select time range and shift it to start from zero;
  • Resample Time, Inputs, Outputs in unique $\delta_T$;
  • Smooth Input/Observe data if flag bit0 is set;
  • Take derivatives of observe data if flag bit1 is set.

In [36]:
# Pick up focused time ranges
time_marks = [[341.47, 350, "free drop 1"],
              [369.83, 377, "free drop 2"],
              [392.38, 400, "free drop 3"],
              [413.57, 421, "free drop 4"],
              [439.80, 450, "free drop 5"],
              [463.46, 474, "free drop 6"],
              [488.48, 498, "free drop 7"],
              [509.40, 520, "free drop 8"],
              [531.80, 543, "free drop 9"],
              [557.84, 569, "free drop 10"],
              [587.60, 599, "free drop 11"],
              [615.00, 627, "free drop 12"]
             ]


# Decide DT,U,Z and their processing method
DT=0.01
process_set = {
            'U':[],
            'Z':[(T_rig, phi_rig,3),],
            'cutoff_freq': 10 #Hz
           }

U_names = []
Y_names = Z_names = ['$\phi_{a,rig} \, / \, ^o$',
                    '$\dot{\phi}_{a,rig} \, / \, ^o/s$',]

display_data_prepare()


$$\delta_T=10.0(ms)$$
Time ranges
IdxStart(s)Spam(s)Description
0341.478.53free drop 1
1369.837.17free drop 2
2392.387.62free drop 3
3413.577.43free drop 4
4439.8010.20free drop 5
5463.4610.54free drop 6
6488.489.52free drop 7
7509.4010.60free drop 8
8531.8011.20free drop 9
9557.8411.16free drop 10
10587.6011.40free drop 11
11615.0012.00free drop 12
$U$
$Z$
$\phi_{a,rig} \, / \, ^o$
$\dot{\phi}_{a,rig} \, / \, ^o/s$

Resample and filter data in sections


In [37]:
resample(True);

Define dynamic model to be estimated

$$\left\{\begin{matrix} \ddot{\phi}_{rig} = \frac{M_{x,rig}}{I_{xx,rig}} \\ M_{x,rig} = M_{x,f} + M_{x,cg} \\ M_{x,f} = -F_c \, sign(\dot{\phi}_{rig}) - f\dot{\phi}_{rig} \\ M_{x,cg} = -m_T g l_{zT} \sin \left ( \phi - \phi_0 \right ) \\ \end{matrix}\right.$$

In [38]:
%%px --local
#update common const parameters in all engines

#problem size
Nx = 2
Nu = 0
Ny = 2
Npar = 7

#reference
S_c = 0.1254   #S_c(m2) 
b_c = 0.7      #b_c(m)
g = 9.81       #g(m/s2)

#other parameters
v_th = 0.5/57.3  #v_th(rad/s)
v_th2 = 0.5/57.3  #v_th(rad/s)

def x0(Z,T,U,params):
    return Z[0,:]/57.3

def xdot(X,t,U,params):
    F_c     = params[0]
    f       = params[1]
    Ixx     = params[2]
    m_T     = params[3]
    l_z_T   = params[4]
    kBrk    = params[5]
    phi0    = params[6]

    phi     = X[0]
    phi_dot = X[1]
        
    moments = -(m_T*l_z_T)*g*math.sin(phi-phi0)

    abs_phi_dot = abs(phi_dot)
    F = f*phi_dot
    if abs_phi_dot > v_th+v_th2:
        F += math.copysign(F_c, phi_dot)
    elif abs_phi_dot > v_th:
        F += math.copysign(F_c*(kBrk-(kBrk-1)*(abs_phi_dot-v_th)/v_th2), phi_dot)
    else:
        F += phi_dot/v_th*(F_c*kBrk)
    moments -= F
    
    phi_dot2 = moments/Ixx
    return [phi_dot, phi_dot2]

def obs(X,T,U,params):
    return X*57.3

In [39]:
display(HTML('<b>Constant Parameters</b>'))
table = ListTable()
table.append(['Name','Value','unit'])
table.append(['$S_c$',S_c,'$m^2$'])
table.append(['$b_c$',b_c,'$m$'])
table.append(['$g$',g,'$m/s^2$'])
display(table)


Constant Parameters
NameValueunit
$S_c$0.1254$m^2$
$b_c$0.7$m$
$g$9.81$m/s^2$

Initial guess

  • Input default values and ranges for parameters
  • Select sections for trainning
  • Adjust parameters based on simulation results
  • Decide start values of parameters for optimization

In [40]:
#initial guess
param0 = [
    0.05,       #F_c(N*m)
    0.03,       #f(N*m/(rad/s))
    0.1678,     #Ixx(kg*m2)
    7.5588,     #m_T(kg)
    0.0526,     #l_z_T(m)
    1.01,       #kBrk
       0,      #phi0(rad)
]

param_name = ['$F_c$','$f$', '$I_{xx,rig}$', '$m_T$',r'$l_{zT}$','$k_{Brk}$','$phi_0$']
param_unit = ['$Nm$',r'$\frac{Nm}{rad/s}$', '$kg\,m^2$', 'kg', 'm','1', 'rad']

NparID = 4
opt_idx = [0,1,2,6]
opt_param0 = [param0[i] for i in opt_idx]
par_del = [0.3*1e-3,  0.3*1e-3,   0.2*1e-3,  0.0174]
bounds =  [(1e-6,0.15),(1e-6,0.1), (1e-6,0.5), (-0.1,0.1)]
display_default_params()

#select sections for training
section_idx = range(8)
display_data_for_train()

#push parameters to engines
push_opt_param()


Default Parameters for optimization
IdxDescription$x_0$unitopt$\delta_{par}$minmax
0$F_c$0.05$Nm$*0.00031e-060.15
1$f$0.03$\frac{Nm}{rad/s}$*0.00031e-060.1
2$I_{xx,rig}$0.1678$kg\,m^2$*0.00021e-060.5
3$m_T$7.5588kg
4$l_{zT}$0.0526m
5$k_{Brk}$1.011
6$phi_0$0rad*0.0174-0.10.1
Training Data
IdxStart(s)Spam(s)Description
0341.478.53free drop 1
1369.837.17free drop 2
2392.387.62free drop 3
3413.577.43free drop 4
4439.8010.20free drop 5
5463.4610.54free drop 6
6488.489.52free drop 7
7509.4010.60free drop 8

In [41]:
# select 2 section from training data
idx = random.sample(section_idx, 2)

interact_guess();


Optimize using ML


In [42]:
display_preopt_params()

if True:
    InfoMat = None
    method = 'trust-ncg'
    def hessian(opt_params, index):
        global InfoMat
        return InfoMat
    dview['enable_infomat']=True
    options={'gtol':1}
    opt_bounds = None
else:
    method = 'L-BFGS-B'
    hessian = None
    dview['enable_infomat']=False
    options={'ftol':1e-4,'maxfun':200}
    opt_bounds = bounds

cnt = 0
tmp_rslt = None
T0 = time.time()
print('#cnt,     Time,      |R|')

%time res =  sp.optimize.minimize(fun=costfunc, x0=opt_param0, \
        args=(opt_idx,), method=method, jac=True, hess=hessian, \
        bounds=opt_bounds, options=options)


Parameters for optimization
IdxDescription$x_0$unitopt$\delta_{par}$minmax
0$F_c$0.05$Nm$*0.00031e-060.15
1$f$0.03$\frac{Nm}{rad/s}$*0.00031e-060.1
2$I_{xx,rig}$0.1678$kg\,m^2$*0.00021e-060.5
3$m_T$7.5588kg
4$l_{zT}$0.0526m
5$k_{Brk}$1.011
6$phi_0$0.0rad*0.0174-0.10.1
#cnt,     Time,      |R|
0001,    0.384,   500256.5
0002,    0.764,   500256.5
0003,    1.083,    42950.1
0004,    1.435,      474.8
0005,    1.837,       94.9
0006,    2.195,       36.8
0007,    2.590,       23.9
0008,    2.961,       23.9
0009,    3.262,       23.9
0010,    3.543,       27.4
0011,    3.839,       31.6
0012,    4.167,       35.1
0013,    4.527,       36.4
0014,    4.833,       36.7
0015,    5.150,       36.8
0016,    5.467,       36.8
0017,    5.768,       36.8
0018,    6.069,       36.8
0019,    6.365,       36.8
0020,    6.665,       36.8
0021,    6.950,       36.8
0022,    7.248,       36.8
0023,    7.655,       36.8
0024,    8.000,       36.8
0025,    8.322,       36.8
0026,    8.647,       36.8
0027,    8.960,       36.8
0028,    9.266,       36.8
0029,    9.603,       36.8
0030,    9.908,       36.8
0031,   10.213,       36.8
0032,   10.503,       36.8
0033,   10.813,       36.8
0034,   11.142,       36.8
0035,   11.428,       36.8
0036,   11.698,       36.8
0037,   12.017,       36.8
0038,   12.313,       36.8
0039,   12.588,       36.8
0040,   12.887,       36.8
0041,   13.233,       36.8
0042,   13.579,       36.8
0043,   13.897,       36.8
Wall time: 13.9 s

Show and test results


In [43]:
display_opt_params()

# show result
idx = range(len(time_marks))
display_data_for_test();

update_guess();


Parameters from optimization
IdxDescription$x$
0$F_c$0.0582804147899
1$f$0.0406548649377
2$I_{xx,rig}$0.198245063396
6$phi_0$-0.00107203709542
Test Data
IdxStart(s)Spam(s)Description
0341.478.53free drop 1
1369.837.17free drop 2
2392.387.62free drop 3
3413.577.43free drop 4
4439.8010.20free drop 5
5463.4610.54free drop 6
6488.489.52free drop 7
7509.4010.60free drop 8
8531.8011.20free drop 9
9557.8411.16free drop 10
10587.6011.40free drop 11
11615.0012.00free drop 12

In [20]:
toggle_inputs()
button_qtconsole()



In [ ]: