网站优化做网站优化,wordpress支持什么格式视频,网站seo优化书籍,外文网站设计目录 说明源代码运行实例workflow_example_1.mworkflow_example_2.mworkflow_example_3.mworkflow_example_4.m 测试1、 结构体兼容性问题2、append的兼容性问题3、修改后的MARRMoT_model.m 说明
MARRMoT是一个新的水文模型比较框架#xff0c;允许不同概念水文模型结构之间的… 目录 说明源代码运行实例workflow_example_1.mworkflow_example_2.mworkflow_example_3.mworkflow_example_4.m 测试1、 结构体兼容性问题2、append的兼容性问题3、修改后的MARRMoT_model.m 说明
MARRMoT是一个新的水文模型比较框架允许不同概念水文模型结构之间的客观比较。该框架为47个独特的模型结构提供了Matlab代码所有模型结构的标准化参数范围以及每个模型的强大数值实现。该框架提供了大量的文档用户手册和几个工作流脚本给予如何使用该框架的例子。包括以下模型
FLEX-TopoIHACRESGR4JTOPMODELSIMHYDVICLASCAMTCMTANKXINANJIANGHYMODSACRAMENTOMODHYDROLOGHBV-96MCRMSMARNAMHYCYMODELGSM-SOCONTECHOPRMSCLASSICIHM19
MARRMoT框架中的模型期望接近它们所基于的原始模型但差异是不可避免的。在许多情况下一个模型存在多个版本但这些版本并不总是很容易区分的。存在一定程度的模型名称相等性其中单个名称用来指代同一基模型的各种不同版本。一个很好的例子是TOPMODEL其中许多变体存在于基于地形指数的最初概念之上。MARRMoT模型往往基于任何给定模型的旧出版物而不是新出版物以接近原作者的“预期”模型但我们的选择是务实的以实现MARRMoT中可用通量和模型结构的更大多样性。每个模型的描述列出了构成该模型MARRMoT版本基础的论文。
MARRMoT框架以模块化的方式设置并以单个的通量文件作为基本的构建块。模型文件通过指定每个模型的内部工作方式来为其定义一类对象而超类文件则定义了所有模型通用的所有过程。图1显示了工具箱结构的示意图概述。
源代码
降雨径流模型的模块化评估
运行实例
workflow_example_1.m
这个示例工作流程包含一个单一模型对单一集水区的应用示例。它包括5个步骤
数据准备模型选择和设置(设置参数)模型求解器设置模型生成和设置模型运行输出可视化
workflow_example_2.m
这个示例工作流程包含一个单一模型对单一集水区的应用示例。它包括5个步骤
数据准备模型选择和设置(默认参数)模型求解器设置模型生成和设置模型运行输出可视化
workflow_example_3.m
这个示例工作流程包含多个模型对单一集水区的应用示例。它包括6个步骤
数据准备模型选择和设置模型求解器设置 对于列表中的每个模型 5. 模型生成和设置 6. 模型运行 7. 输出可视化
workflow_example_4.m
这个示例工作流程包含了一个单一模型对单一集水区的校准应用示例。它包括7个步骤
数据准备模型选择和设置模型求解器设置和时间步进方案校准设置模型校准校准结果评估输出可视化
测试
原始代码根据作者描述在MATLAB版本9.11.0.1873467R2021b上开发的并使用Octave 6.4.0进行了测试。笔者使用R2017a时出现了几个可能是兼容性的问题。
1、 结构体兼容性问题
动态字段或属性名称必须为字符矢量。 将将string转换为char即可使用结构体动态赋值。
2、append的兼容性问题
错误使用 append (line 38) Wrong number of input arguments for obsolete matrix-based syntax. 在MATLAB中append函数是用来将矩阵或数组的元素添加到另一个矩阵或数组的末尾。如果你想将字符串连接到另一个字符串可以使用cat函数或者简单的加号。
3、修改后的MARRMoT_model.m
classdef MARRMoT_model handle
% Superclass for all MARRMoT models% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See https://www.gnu.org/licenses/ for details.properties% attribute to store whether we are running MATLAB or OctaveisOctave % 1 if were on Octave, 0 if MATLAB% static attributes, set for each models in the model definitionnumStores % number of model storesnumFluxes % number of model fluxesnumParams % number of model parametersparRanges % default parameter rangesJacobPattern % pattern of the Jacobian matrix of model store ODEsStoreNames % Names for the storesFluxNames % Names for the fluxesFluxGroups % Grouping of fluxes (useful for water balance and output)StoreSigns % Signs to give to stores (-1 is a deficit store), assumes all 1 if not given% attributes set at the beginning of the simulation% directly by the usertheta % Set of parametersdelta_t % time stepS0 % initial store valuesinput_climate % vector of input climatesolver_opts % options for numerical solving of ODEs% automatically, based on parameter setstore_min % store minimum valuesstore_max % store maximum values% attributes created and updated automatically throughout a% simulationt % current timestepfluxes % vector of all fluxesstores % vector of all storesuhs % unit hydrographs and still-to-flow fluxessolver_data % step-by-step info of solver used and residualsstatus % 0 model created, 1 simulation endedendmethods% This will run as soon as any model object is createdfunction [obj] MARRMoT_model()obj.isOctave exist(OCTAVE_VERSION, builtin)~0;% if running in Octave, load the optim packageif obj.isOctave; pkg load optim; endendfunction [] set.isOctave(obj, value); obj.isOctave value; end% Set methods with checks on inputs for attributes set by the user:function [] set.delta_t(obj, value)if numel(value) 1 || isempty(value)obj.delta_t value;obj.reset();elseerror(delta_t must be a scalar)endendfunction [] set.theta(obj, value)if numel(value) obj.numParams || isempty(value)obj.theta value(:);obj.reset();elseerror([theta must have int2str(obj.numParams) elements])endendfunction [] set.input_climate(obj, value)if isstruct(value)if isfield(value, delta_t)obj.delta_t value.delta_t;elseif isempty(obj.delta_t)error([delta_t is not in input climate struct: ...add it with obj.delta_t delta_t])endif isfield(value, {precip pet temp})P value.precip/obj.delta_t;Ea value.pet/obj.delta_t;T value.temp/obj.delta_t;obj.input_climate [P(:) Ea(:) T(:)];obj.reset();elseerror([Input climate struct must contain fields: ...precip, pet, temp]);endelseif isnumeric(value)if size(value,2) || isempty(value)obj.input_climate value;obj.reset();elseerror([Input climate must have 3 columns: ...precip, pet, temp]);endelseerror([Input climate must either be a struct or ...a numeric array of 3 columns]);endendfunction [] set.S0(obj, value)if numel(value) obj.numStores || isempty(value)obj.S0 value(:);obj.reset();elseerror([S0 must have int2str(obj.numStores) elements])endendfunction [] set.solver_opts(obj, value)% add options to default onesobj.solver_opts obj.add_to_def_opts(value);obj.reset();end% INIT_ runs before each model run to initialise store limits,% auxiliary parameters etc. it calls INIT which is model specificfunction obj init_(obj)% min and max of storesobj.store_min zeros(obj.numStores,1);obj.store_max inf(obj.numStores,1);% empty vectors of fluxes and storest_end size(obj.input_climate, 1);obj.stores zeros(t_end, obj.numStores);obj.fluxes zeros(t_end, obj.numFluxes);% empty struct with the solver dataobj.solver_data.resnorm zeros(t_end,1);obj.solver_data.solver zeros(t_end,1);if(~obj.isOctave); obj.solver_data.solver categorical(obj.solver_data.solver); end;obj.solver_data.iter zeros(t_end,1);% model specific initialisationobj.init();end% RESET is called any time that a user-specified input is changed% (t, delta_t, input_climate, S0, solver_options) and resets any% previous simulation ran on the object.% This is to prevent human error in analysing results.function obj reset(obj)obj.t []; % current timestepobj.fluxes []; % vector of all fluxesobj.stores []; % vector of all storesobj.uhs []; % unit hydrographs and still-to-flow fluxesobj.solver_data []; % step-by-step info of solver used and residualsobj.status 0; % 0 model created, 1 simulation endedend% ODE approximation with Implicit Euler time-stepping schemefunction err ODE_approx_IE(obj, S)S S(:);delta_S obj.model_fun(S);if obj.t 1; Sold obj.S0(:);else; Sold obj.stores(obj.t-1,:);enderr (S - Sold)/obj.delta_t - delta_S;end % SOLVE_STORES solves the stores ODEs function [Snew, resnorm, solver, iter] solve_stores(obj, Sold)solver_opts obj.solver_opts;% This reduces the tolerance to a fraction of the smallest store,% if stores are very small, with 1E-6 as minimum% (if resnorm_tolerance is 0.1 as default)resnorm_tolerance solver_opts.resnorm_tolerance * min(min(abs(Sold)) 1E-5, 1);% create vectors for each of the three solutions (NewtonRaphon,% fsolve and lsqnonlin), this way if all three run it takes the% best at the end and not the last one.Snew_v zeros(3, obj.numStores);resnorm_v Inf(3, 1);iter_v ones(3,1);% first try to solve the ODEs using NewtonRaphsonif(~obj.isOctave) %if MATLAB[tmp_Snew, tmp_fval] ...NewtonRaphson(obj.ODE_approx_IE,...Sold,...solver_opts.NewtonRaphson);else % if Octave[tmp_Snew, tmp_fval] ...NewtonRaphson_octave(obj.ODE_approx_IE,...Sold,...solver_opts.NewtonRaphson);endtmp_resnorm sum(tmp_fval.^2);Snew_v(1,:) tmp_Snew;resnorm_v(1) tmp_resnorm;% if NewtonRaphson doesnt find a good enough solution, run FSOLVEif tmp_resnorm resnorm_tolerance [tmp_Snew,tmp_fval,~,tmp_iter] ...obj.rerunSolver(fsolve, ... tmp_Snew, ... % recent estimatesSold); % storages at previous time steptmp_resnorm sum(tmp_fval.^2);Snew_v(2,:) tmp_Snew;resnorm_v(2) tmp_resnorm;iter_v(2) tmp_iter;% if FSOLVE doesnt find a good enough solution, run LSQNONLINif tmp_resnorm resnorm_tolerance[tmp_Snew,tmp_fval,~,tmp_iter] ...obj.rerunSolver(lsqnonlin, ... tmp_Snew, ... % recent estimatesSold); % storages at previous time steptmp_resnorm sum(tmp_fval.^2);Snew_v(3,:) tmp_Snew;resnorm_v(3) tmp_resnorm;iter_v(3) tmp_iter;endend% get the best solution[resnorm, solver_id] min(resnorm_v);Snew Snew_v(solver_id,:);iter iter_v(solver_id);if(obj.isOctave)solver solver_id;elsesolvers [NewtonRaphson, fsolve, lsqnonlin];solver solvers(solver_id);endend% RERUNSOLVER Restarts a root-finding solver with different % starting pointsfunction [ Snew, fval, stopflag, stopiter ] ...rerunSolver( obj,...solverName,...initGuess,...Sold)% get out useful attributessolver_opts obj.solver_opts.(solverName);solve_fun obj.ODE_approx_IE;max_iter obj.solver_opts.resnorm_maxiter;resnorm_tolerance obj.solver_opts.resnorm_tolerance * min(min(abs(Sold)) 1E-5, 1);% Initialize iteration counter, sampling checker and find number of ODEsiter 1;resnorm resnorm_tolerance 1; % i.e. greater than the required accuracynumStores obj.numStores;stopflag 1; % normal function run% Initialise vector of Snew and fval for each iteration, this way you can% keep the best one, not the last one.Snew_v zeros(numStores, max_iter);fval_v inf(numStores,max_iter);resnorm_v inf(1, max_iter);Snew -1 * ones(numStores, 1);% Start the re-sampling% Re-sampling uses different starting points for the solver to see if% solution accuracy improves. Starting points are alternated as follows:% 1. location where the solver got stuck% 2. storages at previous time step% 3. minimum values% 4. maximum values% 5. randomized values close to solution of previous time stepswhile resnorm resnorm_tolerance% Select the starting pointsswitch itercase 1x0 initGuess(:); % 1. Location where solver got stuckcase 2x0 Sold(:); % 2. Stores at t-1case 3x0 max(-2*10^4.*ones(numStores,1),obj.store_min(:)); % 3. Low values (store minima or -2E4)case 4x0 min(2*10^4.*ones(numStores,1),obj.store_max(:)); % 4. High values (store maxima or 2E4)otherwisex0 max(zeros(numStores,1),...Sold(:)randn(numStores,1).*Sold(:)/10); % 5. Randomized values close to starting locationend% Re-run the solverif strcmpi(solverName, fsolve)[Snew_v(:,iter), fval_v(:,iter), stopflag] ...fsolve(solve_fun, x0, solver_opts);elseif strcmpi(solverName, lsqnonlin)[Snew_v(:,iter), ~, fval_v(:,iter), stopflag] ...lsqnonlin(solve_fun, x0, obj.store_min, [], solver_opts);elseerror(Only fsolve and lsqnonlin are supported);endresnorm_v(iter) sum(fval_v(:,iter).^2);[resnorm,stopiter] min(resnorm_v);fval fval_v(:,stopiter);Snew Snew_v(:,stopiter);% Break out of the loop of iterations exceed the specified maximumif iter max_iterstopflag 0; % function stopped due to iteration countbreakend% Increase the iteration counteriter iter 1;endend% RUN runs the model with a given climate input, initial stores,% parameter set and solver settings.% none of the arguments are needed, they can be set beforehand with% obj.theta theta; obj.input_climate input_climate; etc.% then simply obj.run() without arguments.function [] run(obj,...input_climate,... S0,...theta,...solver_opts)if nargin 4 ~isempty(solver_opts)obj.solver_opts solver_opts;endif nargin 3 ~isempty(theta)obj.theta theta;endif nargin 2 ~isempty(S0)obj.S0 S0;endif nargin 1 ~isempty(input_climate)obj.input_climate input_climate;end% run INIT_ method, this will calculate all auxiliary parameters% and set up routing vectors and store limitsobj.init_();t_end size(obj.input_climate, 1);for t 1:t_endobj.t t;if t 1; Sold obj.S0(:);else; Sold obj.stores(t-1,:);end[Snew,resnorm,solver,iter] obj.solve_stores(Sold);[dS, f] obj.model_fun(Snew);obj.fluxes(t,:) f * obj.delta_t;obj.stores(t,:) Sold dS * obj.delta_t;obj.solver_data.resnorm(t) resnorm;obj.solver_data.solver(t) solver;obj.solver_data.iter(t) iter;obj.step();endobj.status 1;end% GET_OUTPUT runs the model exactly like RUN, but output is% consistent with current MARRMoTfunction [fluxOutput,...fluxInternal,...storeInternal,...waterBalance,...solverSteps] get_output(obj,...varargin)if nargin 1 || isempty(obj.status) || obj.status 0 obj.run(varargin{:});end% --- Fluxes leaving the model --- fg fieldnames(obj.FluxGroups);fluxOutput struct();for k1:numel(fg)idx abs(obj.FluxGroups.(fg{k}));signs sign(obj.FluxGroups.(fg{k}));fluxOutput.(fg{k}) sum(signs.*obj.fluxes(:,idx),2);end% --- Fluxes internal to the model ---fluxInternal struct;for i 1:obj.numFluxesfluxInternal.(char(obj.FluxNames{i})) obj.fluxes(:,i);end% --- Stores ---storeInternal struct;for i 1:obj.numStoresstoreInternal.(char(obj.StoreNames{i})) obj.stores(:,i);end% --- Water balance, if requested ---if nargout 4waterBalance obj.check_waterbalance();end% --- step-by-step data of the solver, if requested ---if nargout 5solverSteps obj.solver_data;endend% CHECK_WATERBALANCE returns the waterbalance% like in MARRMoT1, it will print to screenfunction [out] check_waterbalance(obj, varargin)if nargin 1 || isempty(obj.status) || obj.status 0 obj.run(varargin{:});end% Get variablesP obj.input_climate(:,1);fg fieldnames(obj.FluxGroups);OutFluxes zeros(1,numel(fg));for k1:numel(fg) % cumulative of each flow leaving the modelidx abs(obj.FluxGroups.(fg{k}));signs sign(obj.FluxGroups.(fg{k}));OutFluxes(k) sum(sum(signs.*obj.fluxes(:,idx), 1),2);endif isempty(obj.StoreSigns); obj.StoreSigns repelem(1, obj.numStores); enddS obj.StoreSigns(:) .* (obj.stores(end,:) - obj.S0); % difference of final and initial storage for each storeif isempty(obj.uhs); obj.uhs {}; endR cellfun((uh) sum(uh(2,:)), obj.uhs); % cumulative of each flows still to be routed% calculate water balanceout sum(P) - ... % input from precipitationsum(OutFluxes) - ... % all fluxes leaving the model (some may be entering, but they should have a negative sign)sum(dS) - ... % all differences in storagesum(R); % all flows still being routeddisp([Total P ,num2str(sum(P)), mm.])for k 1:numel(fg)disp([Total ,char(fg(k)), ,...num2str(-OutFluxes(k)), mm.])endfor s 1:obj.numStoresif obj.StoreSigns(s) -1ending (deficit store).;elseending.;enddisp([Delta S,num2str(s), ,...num2str(-dS(s)), mm,ending])endif ~isempty(R)disp([On route ,num2str(-sum(R)), mm.])enddisp(-------------)disp([Water balance , num2str(out), mm.])end% GET_STREAMFLOW only returns the streamflow, runs the model if it% hadnt run already.function Q get_streamflow(obj, varargin)if nargin 1 || isempty(obj.status) || obj.status 0obj.run(varargin{:});endQ sum(obj.fluxes(:,obj.FluxGroups.Q),2);end% CALIBRATE uses the chosen algorithm to find the optimal parameter% set, given model inputs, objective function and observed streamflow.% the function chosen in algorithm should have the same inputs and% outputs as MATLABs fminsearch.function [par_opt,... % optimal parameter setof_cal,... % value of objective function at par_optstopflag,... % flag indicating reason the algorithm stoppedoutput] ... % output, see fminsearch for detailcalibrate(obj,...Q_obs,... % observed streamflowcal_idx,... % timesteps to use for model calibrationoptim_fun,... % function to use for optimisation (must have same structure as fminsearch)par_ini,... % initial parameter estimatesoptim_opts,... % options to optim_funof_name,... % name of objective function to useinverse_flag,... % should the OF be inversed?display,... % should I display information about the calibration?varargin) % additional arguments to the objective functionif isempty(obj.input_climate) || isempty(obj.delta_t) ||...isempty(obj.S0) || isempty(obj.solver_opts)error([input_climate, delta_t, S0 and solver_opts ...attributes must be specified before calling ...calibrate.]);end% if the list of timesteps to use for calibration is empty,% use all steps if isempty(cal_idx)cal_idx 1:length(Q_obs);end% use display by defaultif isempty(display)display true;end% use the data from the start to the last value of cal_idx to% run the simulationif islogical(cal_idx); cal_idx find(cal_idx); endinput_climate_all obj.input_climate;obj.input_climate input_climate_all(1:max(cal_idx),:);Q_obs Q_obs(1:max(cal_idx));% if the initial parameter set isnt set, start from mean% values of parameter rangeif isempty(par_ini)par_ini mean(obj.parRanges,2);end% helper function to calculate fitness given a set of% parametersfunction fitness fitness_fun(par)Q_sim obj.get_streamflow([],[],par);fitness (-1)^inverse_flag*feval(of_name, Q_obs, Q_sim, cal_idx, varargin{:});end% display some useful things for the user to make sure they% used the right settingsif displaydisp(---)disp([Starting calibration of model class(obj) .])disp([Simulation will run for timesteps 1- num2str(max(cal_idx)) .])% this is a bit ugly, but it formats the list of cal_idx in% a pretty and concise waycal_idx sort(cal_idx);i 1;previous cal_idx(i);cal_idx_str num2str(previous); while i numel(cal_idx) i i 1; if cal_idx(i)-previous 1 i find(diff(cal_idx(i:end)) ~ 1, 1) i - 1; if isempty(i); i numel(cal_idx); end previous cal_idx(i); cal_idx_str [cal_idx_str -, num2str(previous)]; else previous cal_idx(i); cal_idx_str [cal_idx_str , , num2str(previous)]; end enddisp([Objective function of_name will be calculated in time steps cal_idx_str .])disp([The optimiser optim_fun will be used to optimise the objective function.])disp([Options passed to the optimiser:])disp(optim_opts)disp(All other options are left to their default,)disp(check the source code of the optimiser to find these default values.)disp(---)end[par_opt,... % optimal parameter set at the end of the optimisationof_cal,... % value of the objective function at par_optstopflag,... % flag indicating reason the algorithm stoppedoutput] ... % output, see fminsearch for detailfeval(optim_fun,... % run the optimisation algorithm chosenfitness_fun,... % function to optimise is the fitness functionpar_ini,... % initial parameter setoptim_opts); % optimiser options% if of_cal was inverted, invert it back before returningof_cal (-1)^inverse_flag * of_cal;% reset the whole input climate as it was before the% calibrationobj.input_climate input_climate_all;end% function to return default solver optionsfunction solver_opts default_solver_opts(obj)solver_opts.resnorm_tolerance 0.1; % Root-finding convergence tolerancesolver_opts.resnorm_maxiter 6; % Maximum number of re-runs used in rerunSolversolver_opts.NewtonRaphson optimset(MaxIter, obj.numStores * 10);% if MATLABif(~obj.isOctave)solver_opts.fsolve optimoptions(fsolve,...Display,none,... % Disable display settingsJacobPattern, obj.JacobPattern);solver_opts.lsqnonlin optimoptions(lsqnonlin,... % lsqnonlin settings for cases where fsolve failsDisplay,none,...JacobPattern,obj.JacobPattern,...MaxFunEvals,1000);else % if OCTAVEsolver_opts.fsolve optimset(Display, off);solver_opts.lsqnonlin optimset(Display, off, MaxFunEvals,1000);endend% function to add new solver opts to the default onesfunction solver_opts add_to_def_opts(obj, opts)def_opts obj.default_solver_opts(); if nargin 1 || isempty(opts)solver_opts def_opts;elsedef_fields fieldnames(def_opts);% for each field in the default options (5 at the moment)for k 1:length(def_fields)field def_fields{k};% if the field is not provided, use the default oneif ~isfield(opts, field) || isempty(opts.(field))solver_opts.(field) def_opts.(field);% if the field is provided, and the dafault is a struct,% add the new values to the default structelseif isstruct(def_opts.(field))solver_opts.(field) optimset(def_opts.(field),opts.(field));% if the field is provided, and the dafault is not a struct,% discard the default and use the new valueelsesolver_opts.(field) opts.(field);endendend endend
end
