www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregusermod/weibul.m
function varargout= weibul(m,x,varargin)
%WEIBUL Weibul growth model
%
% WEIBUL implements the Weibul local model in the Model Browser. This
% file can be used as a template for writing your own user-defined models,
% however this file should not be altered as it will cause the Weibul
% model to be altered and it may no longer work in the Model Browser.
%
% The calling syntax function is VARARGOUT= WEIBUL(U,X,VARARGIN).
% Copyright 2000-2014 The MathWorks, Inc. and Ford Global Technologies, Inc.
% model parameters
b= double(m);
if isnumeric(x)
% evaluation of user-defined model is defined here
% function definition
% make sure the expression is vectorised
y = b(1) - (b(1)-b(2)).*exp(-(b(3).*x).^b(4));
% the weibul function is not defined for x < 0
y(x<0) = NaN;
% this line must not be changed
varargout{1} = y;
else
%---------------------------------------------------------------
% x specifies which internal function to evaluate.
% All internal functions are prefixed by 'i_'.
% The model and model parameters are passed to all internal functions
fInternal = str2func(sprintf('i_%s',lower(x)));
[varargout{1:nargout}]= fInternal(m,b,varargin{:});
%---------------------------------------------------------------
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mandatory functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% n= i_numparams(U,b,varargin);
% n= i_nfactors(U,b,varargin);
% [param,OK]= i_initial(U,b,X,Y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% optional functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% i_jacobian(U,b,X) % analytic jacobian
% i_foptions(U,b,fopts) % optimization parameters
% [LB,UB,A,c,nlcon,optparams]= i_constraints(U,b,X,Y) % define parameter constraints
% g= i_nlconstraints(U,b,X,Y) % evaluate nonlinear constraints
% str= i_char(U,b) % display string for function
% str= i_str_func(U,b) % one line summary of function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% optional local regression functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% names= i_rfnames(U,b) % response feature evaluation
% [rf,dG]= i_rfvals(U,b) % response feature evaluation
% p= i_reconstruct(U,b,Yrf,dG,rfuser) % local model reconstruction
%---------------------------------------------------------------
% internal functions begin here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mandatory functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function n= i_nfactors(U,b) %#ok<INUSD,DEFNU>
%I_NFACTORS number of input factors to model
%
% n = i_nfactors(m,b);
n= 1;
function n= i_numparams(U,b) %#ok<INUSD,DEFNU>
%I_NUMPARAMS number of parameters to be estimated
%
% n = i_numparams(m,b);
n= 4;
function [param,OK]= i_initial(U,b,X,Y) %#ok<INUSL,DEFNU>
%I_INITIAL Initial parameter values
%
% [param,OK]= i_initial(m,b); initial parameters independent of data
% [param,OK]= i_initial(m,b,X,Y); initial parameters dependent on data
% Outputs
% param Initial parameter values
% OK model can be fitted
if nargin>2
% data dependent initial parameter values
OK=1;
Y=double(Y);
X=double(X);
yok= Y>0;
Y= Y(yok);
X= X(yok,:);
if ~isempty(Y)
% Initial estimate of alpha from fitting data
mY=max(Y);
r= mY-min(Y);
% find out if there are a number of points very close to the maximum y
% we are then likely to be very close to the asymptote
f = min(length(find(Y-mY > -r*0.01))+1,5);
alpha= (1 + 10^(-f))*mY;
% Initial estimate of alpha and beta
mY=min(Y);
% find out if there are a number of points very close to the maximum y
% we are then likely to be very close to the asymptote
f= min(length(find(Y-mY < r*0.01))+1,5);
beta= (1 - 10^(-f))*mY;
% Get regression matrix
X=X(:);
Xs= [ones(size(X)) log(X)];
ylin= log(-log((alpha-Y)./(alpha-beta)));
% Find gamma and kappa
p= Xs\ylin;
% parameters from linear regression
delta= p(2) ;
kappa= exp(p(1)/delta);
% Returned parameter initial values
param= [alpha beta kappa delta]';
% Bounds from spec.
param(param<=0)= eps;
else
OK=0;
param=[NaN NaN NaN NaN]';
end
OK= OK && isreal(param);
else
% Initial values of parameters which do not depend on fitting data.
% This case must always be supported
param= [2 1 2 5]';
OK=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% optional functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% i_jacobian % analytic jacobian
% i_foptions % optimization parameters
% i_constraints % define parameter constraints
% i_nlconstraints % evaluate nonlinear constraints
% i_char % display string for function
% i_str_func % one line summary of function
% local regression optional functions
% i_rfnames % response feature evaluation
% i_rfvals % response feature evaluation
% i_reconstruct % local model reconstruction
function [LB,UB,A,c,nlcon,optparams]= i_constraints(U,b,varargin) %#ok<INUSD,DEFNU>
%I_CONSTRAINTS constraint definition for least squares fitting
%
% [LB,UB,A,c,nlcon,optparams]= i_constraints(m,b);
% [LB,UB,A,c,nlcon,optparams]= i_constraints(m,b,X,Y);
% Inputs
% m user-defined model
% b model parameter values
% X,Y fitting data
% Outputs
% LB,UB bounds on parameter values
% A,c linear constraints on parameter values
% nlcon number of nonlinear constraints
% optparams optional optimization parameter
% Lower Bounds, Upper Bounds
% make sure these are column vectors
LB=[eps eps eps eps]';
UB=[1e10 1e10 1e10 1e10]';
% only bounds are used for local regression at present
% don't use bounds with linear constraints if you want to use large scale optimization
% Linear A, Linear b (A*b<c)
% this constraint is -alpha+beta<0 => alpha > beta
A= [-1 1 0 0];
c= 0;
% Number of nonlinear constraints
% you must define 'nlconstraints' if nlcon>0
nlcon= 0;
% Note that fmincon will be used if A and b are not empty, or nlcon>0
% otherwise lsqnonlin is used
% optional parameters for cost function
optparams= [];
function g = i_nlconstraints(m,b,X,Y) %#ok<INUSD,DEFNU>
%I_NLCONSTRAINTS evaluate nonlinear constraints
%
% g = i_nlconstraints(m,b,X,Y);
g = [];
function fopts= i_foptions(U,b,fopts) %#ok<INUSL,DEFNU>
%I_FOPTIONS optimization options
%
% fopts = i_foptions(m,b,fopts)
% fopts is the structure for optimization options. See optimset for more
% details. It is only possible to change the options 'Display',
% 'MaxIter', 'TolFun', 'MaxFunEvals' or 'TolX'.
%
% Always base the output on input fopts.
fopts= optimset(fopts,...
'Display','none');
function J= i_jacobian(U,b,x) %#ok<INUSL,DEFNU>
%I_JACOBIAN analytic jacobian
%
% J = i_jacobian(m,b,X);
x = x(:);
% return empty matrix if not defined
J= zeros(length(x),4);
a=b(1);
beta=b(2);
k=b(3);
d=b(4);
% y = b(1) - (b(1)-b(2)).*exp(-(b(3).*x).^b(4));
ekd= exp(-(k.*x).^d);
j2= (a-beta).*(k.*x).^d.*ekd;
J(:,1)= 1-ekd;
J(:,2)= ekd;
J(:,3)= j2.*d./k;
J(:,4)= j2.*log(k.*x);
function c= i_labels(U,b) %#ok<INUSD,DEFNU>
%I_LABELS user parameter names
%
% c = i_labels(m,b);
c={'\alpha','\beta','\kappa','\delta'};
function str= i_char(U,b) %#ok<INUSD,DEFNU>
%I_CHAR display equation string
%
% str = i_char(m,b);
s= get(U,'symbol');
% this can contain TeX expressions supported by HG text objects
str=sprintf('%.3g - (%.3g-%.3g)*exp(-(%.3g*%s)^{%.3g})',b([1 1 2 3]),detex(s{1}), b(4));
function str= i_str_func(U,b,TeX) %#ok<INUSL,DEFNU>
%I_STR_FUNC short name for model
%
% str = i_str_func(m,b);
s= get(U,'symbol');
if nargin==2 || TeX
s= detex(s);
end
% this can contain TeX expressions supported by HG text objects
lab= labels(U);
str= sprintf('%s - (%s - %s)*exp(-(%s*%s)^{%s})',lab{1},lab{1},lab{2},lab{3},s{1},lab{4});
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% used for local regression only
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rname, default] = i_rfnames(U,b) %#ok<INUSD,DEFNU>
%I_RFNAMES response feature names
%
% [rfnames,InitialRFs] = i_rfnames(m,b);
% Outputs
% rfnames Cell array of nonlinear response feature names
% InitialRFs Initial response features (indices) for user-defined
% model
% response feature names
rname= {'inflex'};
if nargout>1
default = [1 2 5 4]; %{'Alpha','Beta','Inflex','Delta'};
end
function [rf,dG]= i_rfvals(U,b) %#ok<INUSL,DEFNU>
%I_RECONSTRUCT nonlinear reconstruction
%
% p = i_reconstruct(m,b,Yrf,dG,rfuser)
% Inputs
% Yrf response feature values
% dG Jacobian of response features with respect to parameters
% rfuser index to the user-defined response features so you can find
% out which response features are which. rfuser(i) = 0 if the
% rf is a parameter.
%
% If all response features are linear in model parameters then you do not
% need to define 'i_reconstruct'.
% this is an example of how to implement a nonlinear response feature
% response feature definition
K = b(3);
D = b(4);
if D>=1
rf= (1/K)*((D-1)/D)^(1/D);
else
rf= NaN;
end
if nargout>1
% Jacobian of response features with respect to model parameters
if D>=1
dG= [0, 0, -((D-1)/D)^(1/D)/K^2,...
1/K*((D-1)/D)^(1/D)*(-1/D^2*log((D-1)/D)+(1/D-(D-1)/D^2)/(D-1))];
else
dG= [0, 0, 1 1];
end
end
function p= i_reconstruct(U,b,Yrf,dG,rfuser) %#ok<INUSL,DEFNU>
%I_RECONSTRUCT nonlinear reconstruction
%
% p = i_reconstruct(m,b,Yrf,dG,rfuser)
% Inputs
% Yrf response feature values
% dG default Jacobian of response features
% rfuser index to the user-defined response features so you can find
% out which response features are which. rfuser(i) = 0 if the
% rf is a parameter.
%
% If all response features are linear in model parameters then you do not
% need to define 'i_reconstruct'.
% reconstruct linear response features using this line
p= Yrf/dG';
% find which response feature is a nonlinear user-defined response feature
f= rfuser>size(p,2);
if ~any(rfuser==3)
% need to use delta (must be > 1) for reconstruction to work
p(:,4)= max(p(:,4),1+16*eps);
p(:,3)= ((p(:,4)-1)./p(:,4)).^(1./p(:,4))./Yrf(:,f);
end
%---------------------------------------------------------------
% i_evalbuild Simulink evaluation block
%---------------------------------------------------------------
function Blk = i_evalbuild(U, b, sys) %#ok<INUSL,DEFNU>
% evalbuild method adds an evaluation Simulink block to
% a Simulink implementation of a model
Blk= add_block('mbcSLModels/LocalMod/weibulEval',[sys,'/weibulEval']);
% break library link
set_param(Blk,'linkstatus','none');
%---------------------------------------------------------------
% i_slrecon Simulink reconstruction block
%---------------------------------------------------------------
function blk= i_slrecon(U, b, LUM, sname) %#ok<INUSL,DEFNU>
% Adds the appropriate reconstruct block
blk= add_block('mbcSLModels/Reconstruct/weibRecon',...
[sname,'/Reconstruct']);
% break library link
set_param(blk,'linkstatus','none');