diff options
Diffstat (limited to 'estpoly.sci')
| -rw-r--r-- | estpoly.sci | 136 |
1 files changed, 136 insertions, 0 deletions
diff --git a/estpoly.sci b/estpoly.sci new file mode 100644 index 0000000..f1ea48d --- /dev/null +++ b/estpoly.sci @@ -0,0 +1,136 @@ +// Estimates Discrete time estpoly model +// y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t) +// Current version uses random initial guess +// Need to get appropriate guess from OE and noise models + +// Authors: Ashutosh,Harpreet,Inderpreet +// Updated(12-6-16) + +//function [theta_estpoly,opt_err,resid] = estpoly(varargin) +function sys = estpoly(varargin) + + [lhs , rhs] = argn(); + if ( rhs < 2 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "estpoly", rhs); + error(errmsg) + end + + z = varargin(1) + if typeof(z) == 'iddata' then + Ts = z.Ts;unit = z.TimeUnit + z = [z.OutputData z.InputData] + elseif typeof(z) == 'constant' then + Ts = 1;unit = 'seconds' + end + if ((~size(z,2)==2) & (~size(z,1)==2)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "estpoly"); + error(errmsg); + end + + if (~isreal(z)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "estpoly"); + error(errmsg); + end + + n = varargin(2) + if (size(n,"*")<5| size(n,"*")>6) then + errmsg = msprintf(gettext("%s: The order and delay matrix [na nb nc nd nf nk] should be of size [5 6]"), "estpoly"); + error(errmsg); + end + + if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then + errmsg = msprintf(gettext("%s: values of order and delay matrix [nb nc nd nf nk] should be nonnegative integer number "), "estpoly"); + error(errmsg); + end + + na = n(1); nb = n(2); nc = n(3); nd = n(4);nf = n(5); + + if (size(n,"*") == 5) then + nk = 1 + else + nk = n(6); + end + + // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + YDATA = z(:,1); + UDATA = z(:,2); + NDATA = size(UDATA,"*"); + function e = G(p,m) + e = YDATA - _oestpolyfun(UDATA,p,na,nb,nc,nd,nf,nk);//_oestpolyfun(UDATA,p,nd,nc,nf,nb,nk); + endfunction + tempSum = na+nb+nc+nd+nf + p0 = linspace(0.0001,0.001,tempSum)'; + [var,errl] = lsqrsolve(p0,G,size(UDATA,"*")); + //disp(errl) + err = (norm(errl)^2); + //disp(err) + opt_err = err; + resid = G(var,[]); + x = var +// b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff"); +// c = poly([1; var(nb+1:nb+nc)]',"q","coeff"); +// d = poly([1; var(nb+nc+1:nb+nc+nd)]',"q","coeff"); +// f = poly([1; var(nb+nd+nc+1:nd+nc+nf+nb)]',"q","coeff"); + a = poly([1; x(1:na)]',"q","coeff"); + b = poly([repmat(0,nk,1);x(na+1:na+nb)]',"q","coeff"); + c = poly([1; x(na+nb+1:na+nb+nc)]',"q","coeff"); + d = poly([1; x(na+nb+nc+1:na+nb+nc+nd)]',"q","coeff"); + f = poly([1; x(na+nb+nd+nc+1:na+nd+nc+nf+nb)]',"q","coeff"); + t = idpoly(coeff(a),coeff(b),coeff(c),coeff(d),coeff(f),Ts) + + //t = sys;//idpoly(1,coeff(b),coeff(c),coeff(d),coeff(f),Ts) + + // estimating the other parameters + [temp1,temp2,temp3] = predict(z,t) + [temp11,temp22,temp33] = pe(z,t) + //pause + estData = calModelPara(temp1,temp11,na+nb+nc+nd+nf) + //pause + t.Report.Fit.MSE = estData.MSE + t.Report.Fit.FPE = estData.FPE + t.Report.Fit.FitPer = estData.FitPer + t.Report.Fit.AIC = estData.AIC + t.Report.Fit.AICc = estData.AICc + t.Report.Fit.nAIC = estData.nAIC + t.Report.Fit.BIC = estData.BIC + t.TimeUnit = unit + sys = t + + //sys.TimeUnit = unit +endfunction + +function yhat = _oestpolyfun(UDATA,x,na,nb,nc,nd,nf,nk)//(UDATA,x,nd,nc,nf,nb,nk) + x=x(:) + q = poly(0,'q') + tempSum = na+nb+nc+nd+nf + // making polynomials + a = poly([1; x(1:na)]',"q","coeff"); + b = poly([repmat(0,nk,1);x(na+1:na+nb)]',"q","coeff"); + c = poly([1; x(na+nb+1:na+nb+nc)]',"q","coeff"); + d = poly([1; x(na+nb+nc+1:na+nb+nc+nd)]',"q","coeff"); + f = poly([1; x(na+nb+nd+nc+1:na+nd+nc+nf+nb)]',"q","coeff"); + bd = coeff(b*d); cf = coeff(c*f); fc_d = coeff(f*(c-a*d)); + if size(bd,"*") == 1 then + bd = repmat(0,nb+nd+1,1) + end + //pause + maxDelay = max([length(bd) length(cf) length(fc_d)]) + yhat = [YDATA(1:maxDelay)] + for k=maxDelay+1:size(UDATA,"*") + bdadd = 0 + for i = 1:size(bd,"*") + bdadd = bdadd + bd(i)*UDATA(k-i+1) + end + + fc_dadd = 0 + for i = 1:size(fc_d,"*") + fc_dadd = fc_dadd + fc_d(i)*YDATA(k-i+1) + end + cfadd = 0 + for i = 2:size(cf,"*") + cfadd = cfadd + cf(i)*yhat(k-i+1) + end + //pause + yhat = [yhat; [ bdadd + fc_dadd - cfadd ]]; + end +endfunction |