0001 function md=cielocontrol(md)
0002
0003
0004
0005 t1=clock;
0006
0007 analysis='control';
0008
0009
0010
0011 SetParameterDefaults;
0012
0013
0014 m=CreateFEMModel(md,'diagnostic_horiz');
0015
0016
0017 indx=1:6:m.uset.gsize; indx=indx(m.tpart);
0018 indy=2:6:m.uset.gsize; indy=indy(m.tpart);
0019
0020
0021 md.dof=m.uset.fsize;
0022
0023
0024 m.rifts=md.rifts;
0025
0026
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
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;
0036 u_g_obs(2:6:m.uset.gsize)=md.vy_obs(m.part)/md.yts;
0037
0038
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
0048 gradients=cell(numparams,1);
0049 for i=1:numparams,
0050 gradients{i}={};
0051 end
0052 c(1).gradients=gradients; clear gradients;
0053
0054
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
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
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
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
0081 disp(' normalizing directions...');
0082 c(n).gradients=Orth(c(n).gradients,c(n-1).gradients);
0083 disp(' done.');
0084
0085
0086 if md.plot
0087 plot_direction;
0088 end
0089
0090
0091 disp(' optimizing along gradient direction...');
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
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
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
0106 parameters=controlconstrain(parameters,m,m.params,md.control_type);
0107 disp(' done.');
0108
0109
0110 if md.plot
0111 plot_newdistribution;
0112 end
0113
0114
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
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'])