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_Exp051_20150612163239.dat.npz'
def loadData():
# Read and parse raw data
global exp_data
exp_data = np.load(filename)
# Select colums
global T_cmp, da1_cmp, da2_cmp, da3_cmp , da4_cmp
T_cmp = exp_data['data33'][:,0]
da1_cmp = exp_data['data33'][:,3]
da2_cmp = exp_data['data33'][:,5]
da3_cmp = exp_data['data33'][:,7]
da4_cmp = exp_data['data33'][:,9]
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[0].plot(T_cmp,da1_cmp,'r', T_cmp,da2_cmp,'g',
T_cmp,da3_cmp,'b', T_cmp,da4_cmp,'m',
picker=1)
ax[1].plot(T_rig,phi_rig, 'b', picker=2)
ax[0].set_ylabel('$\delta \/ / \/ ^o$')
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()
For each section,
In [8]:
# Decide DT,U,Z and their processing method
process_set1 = {
# Pick up focused time ranges
'time_marks' : [
[10.244684487,71.0176928039,"ramp cmp1 u1"],
[75.590251771,136.676120954,"ramp cmp1 d1"],
[138.085246353,198.267991292,"ramp cmp2 u1"],
[202.947148504,263.947709116,"ramp cmp2 d1"],
[265.320787191,326.086083756,"ramp cmp2 d2"],
[328.651258032,387.065040394,"ramp cmp3 u1"],
[391.213644461,451.562391555,"ramp cmp3 d1"],
[454.129343981,513.129351723,"ramp cmp4 u1"],
[515.902028523,577.039776316,"ramp cmp4 d1"],
[580.1785085,679.5971965,"ramp cmp1/3 d1"],
[683.479931261,783.589410569,"ramp cmp1/3 d2"],
[785.543756555,885.058700497,"ramp cmp1/3 u1"],
[889.116441348,987.763826608,"ramp cmp2/4 d1"],
[992.08250855,1089.55268169,"ramp cmp2/4 u1"],
[1098.5355005,1216.24487648,"ramp cmpall u1"],
[1220.27886043,1340.49709824,"ramp cmpall u2"],
[1343.59303659,1464.78924788,"ramp cmpall d1"],
],
'U':[(T_cmp, da1_cmp,0),
(T_cmp, da2_cmp,0),
(T_cmp, da3_cmp,0),
(T_cmp, da4_cmp,0),],
'U_names' : ['$\delta_{a1,cmp} \, / \, ^o$',
'$\delta_{a2,cmp} \, / \, ^o$',
'$\delta_{a3,cmp} \, / \, ^o$',
'$\delta_{a4,cmp} \, / \, ^o$'],
'Z':[(T_rig, phi_rig,1),],
'Z_names' : ['$\phi_{a,rig} \, / \, ^o$'],
'cutoff_freq': 1, #Hz
'consts' : {'DT':0.1, 'id':0, 'V':30}
}
display_data_set(process_set1)
resample(process_set1, append=False);
In [9]:
%%px --local
#update common const parameters in all engines
angles = range(-40,41,5)
angles[0] -= 1
angles[-1] += 1
del angles[angles.index(0)]
angles_num = len(angles)
#problem size
Nx = 0
Nu = 4
Ny = 1
Npar = 4*angles_num+1
#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)
Clda_cmp = -0.315904095782 #Clda_cmp(1/rad)
Clphi_cmp = -0.0131776778575 #Clphi_cmp(1/rad)
#for short
_m_T_l_z_T_g = -(m_T*l_z_T)*g
def obs(Z,T,U,params,consts):
DT = consts['DT']
ID = consts['id']
V = consts['V']
k1 = np.array(params[0:angles_num])
k2 = np.array(params[angles_num:angles_num*2])
k3 = np.array(params[angles_num*2:angles_num*3])
k4 = np.array(params[angles_num*3:angles_num*4])
phi0 = params[-1]
Clda1 = scipy.interpolate.interp1d(angles, Clda_cmp*0.00436332313*k1*angles,assume_sorted=True)
Clda2 = scipy.interpolate.interp1d(angles, Clda_cmp*0.00436332313*k2*angles,assume_sorted=True)
Clda3 = scipy.interpolate.interp1d(angles, Clda_cmp*0.00436332313*k3*angles,assume_sorted=True)
Clda4 = scipy.interpolate.interp1d(angles, Clda_cmp*0.00436332313*k4*angles,assume_sorted=True)
s = T.size
qbarSb = 0.5*1.225*V*V*S_c*b_c
moments_a = qbarSb*(Clda1(U[:s,0])+Clda2(U[:s,1])
+Clda3(U[:s,2])+Clda4(U[:s,3]))
phi = phi0+np.arcsin(np.clip(-moments_a/_m_T_l_z_T_g, -1, 1))
moments_f = np.copysign(F_c, phi);
moments_p = qbarSb*Clphi_cmp*phi;
phi = phi0+np.arcsin(np.clip(-(moments_a+moments_p+moments_f)/_m_T_l_z_T_g, -1, 1))
return (phi*57.3).reshape((-1,1))
In [11]:
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$'])
table.append(['$C_{l \delta a,cmp}$',Clda_cmp,'$rad^{-1}$'])
table.append(['$C_{l \phi,cmp}$',Clphi_cmp,'$rad^{-1}$'])
display(table)
In [12]:
#initial guess
param0 = [1]*(4*angles_num)+[0]
param_name = ['k_{}_{}'.format(i/angles_num+1, angles[i%angles_num]) for i in range(4*angles_num)] + ['$phi_0$']
param_unit = ['1']*(4*angles_num) + ['$rad$']
NparID = Npar
opt_idx = range(Npar)
opt_param0 = [param0[i] for i in opt_idx]
par_del = [0.001]*(4*angles_num) + [0.0001]
bounds = [(0,1.5)]*(4*angles_num) +[(-0.1, 0.1)]
display_default_params()
#select sections for training
section_idx = range(9)
del section_idx[3]
display_data_for_train()
#push parameters to engines
push_opt_param()
In [13]:
# select 4 section from training data
#idx = random.sample(section_idx, 4)
idx = section_idx[:]
interact_guess();
In [14]:
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 [16]:
display_opt_params()
# show result
idx = range(len(sections))
display_data_for_test();
update_guess();
In [17]:
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])
k2 = np.array(params[angles_num:angles_num*2])
k3 = np.array(params[angles_num*2:angles_num*3])
k4 = np.array(params[angles_num*3:angles_num*4])
Clda_cmp1 = Clda_cmp*0.00436332313*k1*angles
Clda_cmp2 = Clda_cmp*0.00436332313*k2*angles
Clda_cmp3 = Clda_cmp*0.00436332313*k3*angles
Clda_cmp4 = Clda_cmp*0.00436332313*k4*angles
print('angeles = ')
print(angles)
print('Clda_cmpx = ')
print(np.vstack((Clda_cmp1,Clda_cmp2,Clda_cmp3,Clda_cmp4)))
%matplotlib inline
plt.figure(figsize=(12,8),dpi=300)
plt.plot(angles, Clda_cmp1, 'r')
plt.plot(angles, Clda_cmp2, 'g')
plt.plot(angles, Clda_cmp3, 'b')
plt.plot(angles, Clda_cmp4, 'm')
plt.xlabel('$\delta_{a,cmp}$')
plt.ylabel('$C_{l \delta a,cmp}$')
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 [ ]: