cielocontrol

PURPOSE ^

inverse control method solution sequence.

SYNOPSIS ^

function md=cielocontrol(md)

DESCRIPTION ^

    inverse control method solution sequence.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function md=cielocontrol(md)
0002 %    inverse control method solution sequence.
0003 
0004     %timing
0005     t1=clock;
0006     
0007     analysis='control';
0008 
0009 
0010     %configure, and initialise data model and parameter defaults:
0011     SetParameterDefaults;
0012 
0013     %Build finite element model
0014     m=CreateFEMModel(md,'diagnostic_horiz');
0015     
0016     %Build partitioning vectors.
0017     indx=1:6:m.uset.gsize; indx=indx(m.tpart);
0018     indy=2:6:m.uset.gsize; indy=indy(m.tpart);
0019 
0020     % figure out number of dof: just for information purposes.
0021     md.dof=m.uset.fsize; %biggest dof number
0022 
0023     %rifts
0024     m.rifts=md.rifts;
0025 
0026     %create y_s0, for adjoint vector solution.
0027     if ~isempty(m.y_s),
0028         m.y_s0=sparse(size(m.y_s,1),1);
0029     else
0030         m.y_s0=[];
0031     end
0032 
0033     %initialize observations
0034     u_g_obs=zeros(m.uset.gsize,1); 
0035     u_g_obs(1:6:m.uset.gsize)=md.vx_obs(m.part)/md.yts; %from m/a to m/s
0036     u_g_obs(2:6:m.uset.gsize)=md.vy_obs(m.part)/md.yts;
0037 
0038     %initialize parameters on which inversion is done.
0039     numparams=length(md.control_type);
0040     parameters=cell(numparams,1);
0041     for i=1:length(parameters), 
0042         parameter=zeros(m.uset.gsize,1);
0043         parameter(1:6:m.uset.gsize)=eval(['md.' md.control_type{i} '(m.part)']);
0044         parameters{i}=parameter;
0045     end
0046 
0047     %initialize  parameter gradients on which inversion is done.
0048     gradients=cell(numparams,1);
0049     for i=1:numparams,
0050         gradients{i}={};
0051     end
0052     c(1).gradients=gradients; clear gradients;
0053 
0054     %set optimization options.
0055     options=optimset('fminbnd');      
0056     options=optimset(options,'Display','iter');
0057     options=optimset(options,'MaxIter',100);
0058     options=optimset(options,'TolX',md.tolx);
0059     options=optimset(options,'FunValCheck','on');
0060 
0061     for n=2:md.nsteps+1,
0062     
0063         %set options
0064         options=optimset(options,'MaxFunEvals',md.maxiter(n-1));
0065 
0066         disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n-1) '/' num2str(md.nsteps)]));
0067 
0068         %initialize inputs, ie parameters on which we invert.
0069         for i=1:numparams,
0070             eval(['inputs.' md.control_type{i} '=parameters{' num2str(i) '};']);
0071         end
0072         inputs.fit=md.fit(n-1);
0073 
0074         %call gradJ module
0075         disp('      computing gradJ...');
0076         c(n).gradients=GradJCompute(m,m.params,inputs,u_g_obs,md.control_type, analysis);
0077         disp('      done.');
0078 
0079         
0080         %normalize directions.
0081         disp('      normalizing directions...');
0082         c(n).gradients=Orth(c(n).gradients,c(n-1).gradients);
0083         disp('      done.');
0084 
0085         %visualize direction.
0086         if md.plot
0087             plot_direction;
0088         end
0089 
0090         %search along the direction for the B vector that minimizes the misfit J
0091         disp('      optimizing along gradient direction...'); %we use fmincon, see help fmincon.
0092         if numparams==1,
0093             [search_vector c(n).J]=fminbnd('objectivefunctionC',-1,1,options,m,u_g_obs,parameters,md.optscal(n-1,:),md.fit(n-1),c(n).gradients,md.control_type,analysis);
0094         else %optimization of several parameters
0095             [search_vector c(n).J]=fmincon('objectivefunctionC',[0;0],[],[],[],[],[-1;-1],[1;1],[],options,m,u_g_obs,parameters,md.optscal(n-1,:),md.fit(n-1),c(n).gradients,md.control_type,analysis);
0096         end
0097         disp('      done.');
0098     
0099         %update parameters using optimized search vector
0100         for i=1:numparams,
0101             parameters{i}=parameters{i}+search_vector(i)*md.optscal(n-1,i)*c(n).gradients{i};
0102         end
0103 
0104         disp('      constraining the new distribution...');    
0105         %constrain control
0106         parameters=controlconstrain(parameters,m,m.params,md.control_type);
0107         disp('      done.');
0108     
0109         %visualize direction.
0110         if md.plot
0111             plot_newdistribution;
0112         end
0113 
0114         %some temporary saving
0115         if(mod(n,5)==0),
0116             solution=controlfinalsol(c,m,parameters,md.control_type,analysis);
0117             save temporary_control_results solution
0118         end
0119         
0120         disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
0121 
0122     end
0123 
0124     %Create final solution
0125     disp('      preparing final velocity solution...');
0126     solution=controlfinalsol(c,m,parameters,md.control_type,analysis);
0127     disp('      done.');
0128     
0129     disp('      load results ontol model...');
0130     md=loadcontrolfinalsol(md,solution);
0131     disp('      done.');
0132     
0133     
0134     t2=clock;
0135     disp(['      overall time spent on control code ' num2str(etime(t2,t1)) ' seconds'])

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003