%--------------------------------------------------------------------------------------------------% % Copyright (C) <2011> by % % % Permission is hereby granted, free of charge, to any person obtaining a copy of this software and % associated documentation files (the "Software"), to deal in the Software without restriction, % including without limitation the rights to use, copy, modify, merge, publish, distribute, % sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is % furnished to do so, subject to the following conditions: % % The above copyright notice and this permission notice shall be included in all copies or % substantial portions of the Software. % % THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT % NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND % NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, % DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT % OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. %--------------------------------------------------------------------------------------------------% %--------------------------------------------------------------------------------------------------% % % main program for optimization % main_opt.m calls the following functions/scripts in sequence: % (1) opt_input, (2) find_taugamma_strain, and (3)get_residual, (4) reconstruct % get_residual calls the hem kernel: hem90_pol_strain %--------------------------------------------------------------------------------------------------% clear all; tic format short e; addpath('utility'); %--------------------------------------------------------------------------------------------------% % Read in user supplied parameters % [opt,ratedata,params,scales,non_opt_params] = opt_input; %global filename dataname opt.csnum; % %--------------------------------------------------------------------------------------------------% % Get experimental data % tmp=load(opt.filename,opt.dataname); input = getfield(tmp,opt.dataname); % if opt.runflag ==2 nstep = 10; % use more points in plotting. else nstep = 20; end data = []; % Define a struct data that contains time, field, polarizaion and strain. % Properly center the experimental data helps the optimization. % Since the polarization and field are not symmetric, we recenter them so % that they are rotational symmetric. diffp = max(input.P)+min(input.P); %shift in polarization if(opt.csnum==1) diffe = 0.525-0.47; %shift in input field data.T = input.T(1901:nstep:end); %time data.E = (input.E(1901:nstep:end)+0.5*diffe)*1e6; % input field data.Sigma = non_opt_params.SIGMA; %prestress data.P = input.P(1901:nstep:end)-0.5*diffp; %polarization data data.S = input.S(1901:nstep:end)/100; %displacement data else diffe = 0.5154-0.4642; %shift in input field data.T = input.T(1:nstep:end); data.E = (input.E(1:nstep:end)+0.5*diffe)*1e6; data.Sigma = non_opt_params.SIGMA; data.P = input.P(1:nstep:end)-0.5*diffp; %polarization data data.S = input.S(1:nstep:end)/100; %displacement data end % %--------------------------------------------------------------------------------------------------% %Initialize parameters if opt.runflag == 2 % load final estimates to plot if(opt.csnum==1) load results/estimated_values_zerostress; else load results/estimated_values_prestress; end params_opt0 = [params.P_R_p,params.P_R_90,params.EPS_R_p,params.EPS_R_90,params.gamma,params.tau_90,... params.chi_p,params.d_p,params.d_90,params.s_p,... params.alpha_int,params.alpha_frac,params.beta_int,params.beta_frac]; scales_opt0 = [scales.P_R_p,scales.P_R_90,scales.EPS_R_p,scales.EPS_R_90,scales.gamma,scales.tau_90,... scales.chi_p,scales.d_p,scales.d_90,scales.s_p,... scales.alpha_int,scales.alpha_frac,scales.beta_int,scales.beta_frac]; sprintf('The estimated parameters are') params_opt0.*scales_opt0 else % % Initialization % % Call subroutine find_taugamma_strain.m to identify tau and gamma % % Defind rate-dependent data lb1 = ratedata.lb1; ub1 = ratedata.ub1; lb2 = ratedata.lb2; ub2 = ratedata.ub2; step1 = ratedata.step1; step2 = ratedata.step2; ratedata1.T =input.T(lb1:nstep:ub1); ratedata1.E =(input.E(lb1:nstep:ub1)+0.5*diffe)*1e6; ratedata1.Sigma =non_opt_params.SIGMA; ratedata1.P =input.P(lb1:nstep:ub1)-0.5*diffp; ratedata1.S =input.S(lb1:nstep:ub1)/100; ratedata1.step = step1; ratedata2.T =input.T(lb2:nstep:ub2); ratedata2.E =(input.E(lb2:nstep:ub2)+0.5*diffe)*1e6; ratedata2.Sigma =non_opt_params.SIGMA; ratedata2.P =input.P(lb2:nstep:ub2)-0.5*diffp; ratedata2.S =input.S(lb2:nstep:ub2)/100; ratedata2.step = step2; % Use strain to identify tau and gamma [gamma,tau_90]=find_taugamma_strain(ratedata1,ratedata2,params,scales,non_opt_params); % Rescale gamma and tau params.gamma = gamma/scales.gamma; params.tau_90 = tau_90/scales.tau_90; % Initial parameters without the scaling. params_opt0 = [params.P_R_p,params.P_R_90,params.EPS_R_p,params.EPS_R_90,params.gamma,params.tau_90,... params.chi_p,params.d_p,params.d_90,params.s_p,... params.alpha_int,params.alpha_frac,params.beta_int,params.beta_frac]; % Scaling vector scales_opt0 = [scales.P_R_p,scales.P_R_90,scales.EPS_R_p,scales.EPS_R_90,scales.gamma,scales.tau_90,... scales.chi_p,scales.d_p,scales.d_90,scales.s_p,... scales.alpha_int,scales.alpha_frac,scales.beta_int,scales.beta_frac]; sprintf('The initial parameters are') params_opt0.*scales_opt0 end % %--------------------------------------------------------------------------------------------------% % if opt.runflag == 1 save('results/initial_values.mat','params','scales','non_opt_params'); end % %--------------------------------------------------------------------------------------------------% % Call optimization scheme % The residual function f = @(params_opt)get_residual(params_opt,scales_opt0,non_opt_params,data,opt); % lower bound of parameters lb = zeros(1,length(params_opt0)); % options to run optimization without analytic Jacobian % options = optimset('DiffMinChange',1e-5,'Display','iter','TolFun',1e-4, ... % 'Algorithm','trust-region-reflective','ScaleProblem','Jacobian',... % 'TypicalX',scales_opt0); % options to run optimization with analytic Jacobian options = optimset('Display','iter','TolFun',1e-4, ... 'Algorithm','trust-region-reflective','Jacobian','on'); [params_est,resnorm,resid]=lsqnonlin(f,params_opt0,lb,[],options); % The first point used in calculating the residual nn = opt.int; l2_relatvie = sqrt(sum(resid.^2)/(sum(data.P(nn:end).^2)+sum((data.S(nn:end)*opt.wt).^2))) %--------------------------------------------------------------------------------------------------% % save the estimated parameters. params=reconstruct(params_est,params,non_opt_params); if opt.runflag ==0 if(opt.csnum==1) save('results/estimated_values_zerostress.mat','params','scales','non_opt_params'); else save('results/estimated_values_prestress.mat','params','scales','non_opt_params'); end sprintf('The estimated parameters are') params_est.*scales_opt0 end toc