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
In [3]:
filename = 'FIWT_Exp036_20150605144438.dat.npz'
def loadData():
# Read and parse raw data
global exp_data
exp_data = np.load(filename)
# Select colums
global T_cmp, da_cmp
T_cmp = exp_data['data33'][:,0]
da_cmp = exp_data['data33'][:,3]
global T_rig, phi_rig
T_rig = exp_data['data44'][:,0]
phi_rig = exp_data['data44'][:,2]
loadData()
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 [4]:
def checkInputOutputData():
#check inputs/outputs
fig, ax = plt.subplots(2,1,True)
ax[1].plot(T_cmp,da1_cmp,'r', picker=2)
ax[0].plot(T_rig,phi_rig, 'b', picker=2)
ax[1].set_ylabel('$\delta \/ / \/ ^o$')
ax[0].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()
In [16]:
process_set1 = {
'time_marks' : [
[17.6940178696,117,"ramp u1"],
[118.7,230.395312673,"ramp d1"],
[258.807481992,357.486188688,"ramp u2"],
[359.463122988,459.817499014,"ramp d2"],
[461.067939262,558.784538108,"ramp d3"],
[555.553175853,658.648739191,"ramp u4"],
],
'Z':[(T_cmp, da_cmp,1),],
'Z_names' : ['$\delta_{a,cmp} \, / \, ^o$'],
'U_names' : ['$\phi_{a,rig} \, / \, ^o$'],
'U':[(T_rig, phi_rig,1),],
'cutoff_freq': 1, #Hz
'consts':{'DT':0.5,'id':0,'V':30}
}
display_data_set(process_set1)
resample(process_set1, append=False);
In [17]:
%%px --local
#update common const parameters in all engines
angles = range(-170,171,10)
angles_num = len(angles)
#problem size
Nx = 0
Nu = 1
Ny = 1
Npar = angles_num
#reference
S_c = 0.1254 #S_c(m2)
b_c = 0.7 #b_c(m)
g = 9.81 #g(m/s2)
#static measurement
m_T = 7.5588 #m_T(kg)
l_z_T = 0.0424250531303 #l_z_T(m)
#previous estimations
F_c = 0.0532285873599 #F_c(N*m)
#for short
_m_T_l_z_T_g = -(m_T*l_z_T)*g
angles_cmpx = [-41, -35, -30, -25, -20, -15, -10, -5, 5, 10, 15, 20, 25, 30, 35, 41]
Clda_cmpx = np.array([[ 0.05968291, 0.0553231, 0.04696643, 0.03889542, 0.03162213, 0.0271057,
0.01839293, 0.00900769, -0.00497219, -0.01407865, -0.02279645, -0.02607804,
-0.03158932, -0.03814162, -0.04606267, -0.05423365],
[ 0.04821761, 0.0404353, 0.03394095, 0.0276883, 0.0223423, 0.01843785,
0.01272706, 0.00739553, -0.00449509, -0.00899431, -0.01628838, -0.01994654,
-0.02522128, -0.03103204, -0.03819486, -0.04861837],
[ 0.04235569, 0.03650224, 0.02990824, 0.02505072, 0.01880923, 0.01535756,
0.01121967, 0.00679852, -0.0041091, -0.00866967, -0.01516793, -0.01995755,
-0.02558169, -0.03024035, -0.0346464, -0.03949022],
[ 0.03393535, 0.03038317, 0.02596484, 0.02259024, 0.01868117, 0.01450576,
0.00870284, 0.00571733, -0.00445593, -0.00858266, -0.01293602, -0.0150844,
-0.01936097, -0.02433856, -0.03052023, -0.0351404 ]])
Clda_cmp = np.sum(Clda_cmpx, axis=0)
Da_cmp = scipy.interpolate.interp1d(Clda_cmp, angles_cmpx)
def obs(Z,T,U,params,consts):
DT = consts['DT']
ID = consts['id']
V = consts['V']
unk_moment = scipy.interpolate.interp1d(angles, params[0:angles_num])
s = T.size
phi = U[:s,0]
hdr = int(1/DT)
col_fric = np.copysign(-F_c,phi[hdr:]-phi[:-hdr])
col_fric = np.concatenate((np.ones(hdr-1)*col_fric[0], col_fric, np.ones(hdr-1)*col_fric[-1]))
qbarSb = 0.5*1.225*V*V*S_c*b_c
Da = Da_cmp((-_m_T_l_z_T_g*np.sin(phi/57.3)-unk_moment(phi)+col_fric)/qbarSb)
return Da.reshape((-1,1))
In [18]:
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$'])
table.append(['$m_T$',m_T,'$kg$'])
table.append(['$l_{zT}$',l_z_T,'$m$'])
table.append(['$F_c$',F_c,'$Nm$'])
display(table)
In [19]:
#initial guess
param0 = [0]*angles_num
param_name = ['Mu_{}'.format(angles[i]) for i in range(angles_num)]
param_unit = ['Nm']*angles_num
NparID = Npar
opt_idx = range(Npar)
opt_param0 = [param0[i] for i in opt_idx]
par_del = [0.01]*(angles_num)
bounds = [(-2,2)]*(angles_num)
display_default_params()
#select sections for training
section_idx = range(len(sections))
display_data_for_train()
#push parameters to engines
push_opt_param()
In [20]:
# select 4 section from training data
#idx = random.sample(section_idx, 4)
idx = section_idx[:]
# interact_guess();
%matplotlib inline
update_guess();
In [21]:
display_preopt_params()
if False:
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-6,'maxfun':400}
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)
In [23]:
display_opt_params()
# show result
idx = range(len(sections))
display_data_for_test();
update_guess();
In [24]:
res_params = res['x']
params = param0[:]
for i,j in enumerate(opt_idx):
params[j] = res_params[i]
k1 = np.array(params[0:angles_num])
print('angeles = ')
print(angles)
print('L_unk = ')
print(k1)
%matplotlib inline
plt.figure(figsize=(12,8),dpi=300)
plt.plot(angles, k1, 'r')
plt.ylabel('$L_{unk} \, / \, Nm$')
plt.xlabel('$\phi \, / \, ^o/s$')
plt.show()
In [14]:
toggle_inputs()
button_qtconsole()
In [15]:
(-0.05-0.05)/(80/57.3)
Out[15]:
In [16]:
Clda_cmp/4
Out[16]:
In [ ]: