Commit c2a37b2b authored by Nicolas Rodriguez's avatar Nicolas Rodriguez

Upload New File

parent c6bcd1b8
%
% Copyright 2019-2020 LIST (Luxembourg Institute of Science and Technology), all right reserved.
%
% Author: Nicolas Rodriguez (nicolas.bjorn.rodriguez@gmail.com)
%
function [Eff_C,model_output,model_output_warmup,options]=model_run_virtual(param, model_input, meas, options)
% Warm-up (here, it is not used)
if options.nbyrs_warmup>0
if options.display==1
disp('Spinning-up the model...');
tic
end
% Warm-up
options.periods.calc=options.periods.warmup;
y0=options.y0;
[model_output_warmup,options]=model_calc(model_input,param,y0, options);
if model_output_warmup.feas==1
options.ages.init_RTD=model_output_warmup.ages.store_ps(:,end)';
y0=[model_output_warmup.S(end);model_output_warmup.C(end)];
else
model_output.feas=0;
model_output.errtypes=model_output_warmup.errtypes;
Eff_C=-Inf;
return;
end
if options.display==1
toc()
end
else
model_output_warmup.feas=1;
model_output_warmup.C=[];
model_output_warmup.C_S=[];
model_output_warmup.C_Q=[];
model_output_warmup.C_ET=[];
model_output_warmup.S=[];
model_output_warmup.ET=[];
model_output_warmup.M=[];
end
if options.display==1
disp('Running the model...');
tic
end
% Run
if model_output_warmup.feas==1
if exist('y0','var')==0
y0=options.y0;
end
options.periods.calc=options.periods.run;
[model_output,options]=model_calc(model_input,param, y0,options);
if options.display==1
toc()
end
if model_output.feas==0 % if there are numerical issues, the efficiency is -Inf
Eff_C=-Inf;
else
% Assignment of the measured variables
t=model_input.t;
Q=model_input.Q;
Cq=meas.Cq;
% Assignment of the output variables
Csim=cat(2,model_output_warmup.C_Q,model_output.C_Q(options.periods.run));
C=cat(2,model_output_warmup.C,model_output.C(options.periods.run));
C_S=cat(2,model_output_warmup.C_S,model_output.C_S(options.periods.run));
% Calculation of a customized likelihood
period_eval=options.periods.eval;
% For tracer, interpolate simulated values to the timestamps of the observed ones
x1=find(Cq(:,1)>=t(period_eval(1)));
x2=find(Cq(x1,1)<=t(period_eval(end)));
x=x1(x2);
if ~isempty(x)
obs=Cq(x,2);
sim = interp1(t(period_eval),Csim(period_eval),Cq(x,1));
else
obs=[];
sim=[];
end
% evaluate performance
[NS,~]=NSE_virtual(obs,sim);
Eff_C=NS;
end
end
end
\ No newline at end of file
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment