diff options
| -rw-r--r-- | ar.sci | 61 | ||||
| -rw-r--r-- | armaX.sci | 48 | ||||
| -rw-r--r-- | arx.sci | 75 | ||||
| -rw-r--r-- | bj.sci | 64 | ||||
| -rw-r--r-- | buildmacros.sce | 10 | ||||
| -rw-r--r-- | calModelPara.sci | 2 | ||||
| -rw-r--r-- | cleanmacros.sce | 20 | ||||
| -rw-r--r-- | compare.sci | 250 | ||||
| -rw-r--r-- | dataSlice.sci | 54 | ||||
| -rw-r--r-- | estpoly.sci | 89 | ||||
| -rw-r--r-- | etfe.sci | 44 | ||||
| -rw-r--r-- | fitValue.sci | 14 | ||||
| -rw-r--r-- | fitch.sci | 44 | ||||
| -rw-r--r-- | frd.sci | 43 | ||||
| -rw-r--r-- | frdplot.sci | 19 | ||||
| -rwxr-xr-x | iddata.sci | 31 | ||||
| -rw-r--r-- | iddataplot.sci | 44 | ||||
| -rw-r--r-- | identTime.sci | 23 | ||||
| -rw-r--r-- | idinput.sci | 213 | ||||
| -rwxr-xr-x | idpoly.sci | 58 | ||||
| -rw-r--r-- | impulseest.sci | 35 | ||||
| -rw-r--r-- | iv.sci | 61 | ||||
| -rw-r--r-- | iv4.sci | 71 | ||||
| -rw-r--r-- | misdata.sci | 36 | ||||
| -rw-r--r-- | nInputData.sci | 43 | ||||
| -rw-r--r-- | oe.sci | 68 | ||||
| -rw-r--r-- | pe.sci | 50 | ||||
| -rw-r--r-- | predict.sci | 75 | ||||
| -rw-r--r-- | rarx.sci | 40 | ||||
| -rw-r--r-- | resid.sci | 86 | ||||
| -rw-r--r-- | scilab_error.sci | 20 | ||||
| -rw-r--r-- | scilab_sum.sci | 10 | ||||
| -rw-r--r-- | sim.sci | 30 | ||||
| -rw-r--r-- | spa.sci | 30 | ||||
| -rw-r--r-- | stepest.sci | 29 |
35 files changed, 1523 insertions, 367 deletions
@@ -1,11 +1,46 @@ -// Estimates Discrete time AR model -// A(q)y(t) = e(t) -// Current version uses random initial guess -// Authors: Ashutosh,Harpreet,Inderpreet -// Updated(12-6-16) -function sys = ar(varargin) -// +function sys = ar(varargin) +// Parameters Estimation of AR model using Input Output time-domain data +// +// Calling Sequence +// sys = ar(ioData,[na]) +// +// Parameters +// ioData : iddata or [outputData inputData] ,matrix of nx2 dimensions, type plant data +// na : non-negative integer number specified as order of the polynomial A(z^-1) +// sys : idpoly type polynomial have estimated coefficients of A(z^-1) polynomials +// +// Description +// Fit AR model on given input output data +// The mathematical equation of the AR model +// <latex> +// begin{eqnarray} +// A(q)y(t) = e(t) +// end{eqnarray} +// </latex> +// It is SISO type model. It minimizes the sum of the squares of nonlinear functions using Levenberg-Marquardt algorithm. +// sys ,an idpoly type class, have different fields that contains estimated coefficients, sampling time, time unit and other estimated data in Report object. +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.5];b = [0 2 3]; +// model = idpoly(a,b,'Ts',0.1); +// y = sim(u,model) + rand(length(u),1); +// plantData = iddata(y,[],0.1); +// sys = ar(plantData,[2]) +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]); +// a = [1 0.5];b = [0 0.2 0.3]; +// model = idpoly(a,b,'Ts',0.1); +// y = sim(u,model) + rand(length(u),1); +// plantData = [y]; +// sys = ar(plantData,[2]) +// +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + + [lhs , rhs] = argn(); if ( rhs < 2 ) then errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "ar", rhs); @@ -41,7 +76,7 @@ function sys = ar(varargin) end na = n; nb = 0; nk = 0; - // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively YDATA = z(:,1); UDATA = zeros(size(z,1),1) NDATA = size(UDATA,"*"); @@ -59,12 +94,12 @@ function sys = ar(varargin) a = (poly([1,-coeff(a)],'q','coeff')) t = idpoly(coeff(a),1,1,1,1,Ts) - // estimating the other parameters + // estimating the other parameters [temp1,temp2,temp3] = predict([YDATA UDATA],t) [temp11,temp22,temp33] = pe([YDATA UDATA],t) estData = calModelPara(temp1,temp1,n(1)) - //pause + // pause t.Report.Fit.MSE = estData.MSE t.Report.Fit.FPE = estData.FPE t.Report.Fit.FitPer = estData.FitPer @@ -74,15 +109,15 @@ function sys = ar(varargin) t.Report.Fit.BIC = estData.BIC t.TimeUnit = unit sys = t - //sys = idpoly(coeff(a),1,1,1,1,Ts) -// sys.TimeUnit = unit + // sys = idpoly(coeff(a),1,1,1,1,Ts) +// sys.TimeUnit = unit endfunction function yhat = _objfun(UDATA,YDATA,x,na,nb,nk) x=x(:) q = poly(0,'q') tempSum = nb+na - // making polynomials + // making polynomials b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff"); a = 1 - poly([x(nb+1:nb+na)]',"q","coeff") aSize = coeff(a);bSize = coeff(b) @@ -1,5 +1,49 @@ function sys = armaX(varargin) +// Parameters Estimation of ARMAX model using Input Output time-domain data +// +// Calling Sequence +// sys = armaX(ioData,[na nb nc nk]) +// +// Parameters +// ioData : iddata or [outputData inputData] ,matrix of nx2 dimensions, type plant data +// na : non-negative integer number specified as order of the polynomial A(z^-1) +// nb : non-negative integer number specified as order of the polynomial B(z^-1)+1 +// nc : non-negative integer number specified as order of the polynomial C(z^-1) +// nk : non-negative integer number specified as input output delay, Default value is 1 +// sys : idpoly type polynomial have estimated coefficients of A(z^-1),B(z^-1) and C(z^-1) polynomials +// +// Description +// Fit ARMAX model on given input output data +// The mathematical equation of the ARMAX model +// <latex> +// begin{eqnarray} +// A(q)y(n) = B(q)u(n) + C(q)e(n) +// end{eqnarray} +// </latex> +// It is SISO type model. It minimizes the sum of the squares of nonlinear functions using Levenberg-Marquardt algorithm. +// +// sys ,an idpoly type class, have different fields that contains estimated coefficients, sampling time, time unit and other estimated data in Report object. +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.5];b = [0 2 3]; +// model = iddata(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// plantData = iddata(y,u,0.1) +// sys = armaX(plantData,[2,2,1]) +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.5];b = [0 2 3]; +// model = iddata(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// plantData = [y,u] +// sys = armaX(plantData,[2,2,1]) +// +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + [lhs , rhs] = argn(); if ( rhs < 2 ) then errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "armaX", rhs); @@ -34,7 +78,7 @@ function sys = armaX(varargin) error(errmsg); end - na = n(1); nb = n(2); nc = n(3); //nd = n(4);nf = n(5); + na = n(1); nb = n(2); nc = n(3); // nd = n(4);nf = n(5); if (size(n,"*") == 3) then nk = 1 @@ -42,7 +86,7 @@ function sys = armaX(varargin) nk = n(4); end - // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively YDATA = z(:,1); UDATA = z(:,2); @@ -1,28 +1,47 @@ - -// Estimates Discrete time ARX model -// A(q)y(t) = B(q)u(t) + e(t) -// Current version uses random initial guess +function sys = arx(varargin) +// Parameters Estimation of ARX model using Input Output time-domain data +// +// Calling Sequence +// sys = arx(ioData,[na nb nk]) +// +// Parameters +// ioData : iddata or [outputData inputData] ,matrix of nx2 dimensions, type plant data +// na : non-negative integer number specified as order of the polynomial A(z^-1) +// nb : non-negative integer number specified as order of the polynomial B(z^-1)+1 +// nk : non-negative integer number specified as input output delay, Default value is 1 +// sys : idpoly type polynomial have estimated coefficients of A(z^-1) and B(z^-1) polynomials +// +// Description +// Fit ARX model on given input output data +// The mathematical equation of the ARX model +// <latex> +// begin{eqnarray} +// A(q)y(n) = B(q)u(n) + e(t) +// end{eqnarray} +// </latex> +// It is SISO type model. It minimizes the sum of the squares of nonlinear functions using Levenberg-Marquardt algorithm. +// +// sys ,idpoly type, have different fields that contains estimated coefficients, sampling time, time unit and other estimated data in Report object. // - -// Authors: Ashutosh,Harpreet,Inderpreet -// Updated(12-6-16) - // Examples -//loadmatfile("data.mat") -//sys = arx(data,[2,2,1]) -//sys = -// -// A(z) = 1 - 1.3469229 z^-1 + 0.7420890 z^-2 -// -// B(z) = 1.3300766 z^-1 - 0.5726208 z^-2 -// -// Sampling Time = 1 seconds -// -// MSE FPE FitPer AIC AICc nAIC BIC -// 7.4091 7.4388 49.9726 9689.1801 194.2693 2.0067 9711.5838 - +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.5];b = [0 2 3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// plantData = iddata(y,u,0.1) +// sys = arx(plantData,[2,2,1]) +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.5];b = [0 2 3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// plantData = [y,u] +// sys = arx(plantData,[2,2,1]) +// +// Authors +// Ashutosh Kumar Bhargava, Harpreet, Inderpreet -function sys = arx(varargin) [lhs , rhs] = argn(); if ( rhs < 2 ) then errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "arx", rhs); @@ -57,15 +76,15 @@ function sys = arx(varargin) error(errmsg); end - na = n(1); nb = n(2); //nk = n(3); //nf = n(4); -// + na = n(1); nb = n(2); // nk = n(3); // nf = n(4); +// if (size(n,"*") == 2) then nk = 1 else nk = n(3); end - // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively YDATA = z(:,1); UDATA = z(:,2); NDATA = size(UDATA,"*"); @@ -84,12 +103,12 @@ function sys = arx(varargin) a = (poly([1,-coeff(a)],'q','coeff')) t = idpoly(coeff(a),coeff(b),1,1,1,Ts) - // estimating the other parameters + // estimating the other parameters [temp1,temp2,temp3] = predict(z,t) [temp11,temp22,temp33] = pe(z,t) estData = calModelPara(temp1,temp1,n(1)+n(2)) - //pause + // pause t.Report.Fit.MSE = estData.MSE t.Report.Fit.FPE = estData.FPE t.Report.Fit.FitPer = estData.FitPer @@ -106,7 +125,7 @@ function yhat = _objfunarx(UDATA,YDATA,x,na,nb,nk) x=x(:) q = poly(0,'q') tempSum = nb+na - // making polynomials + // making polynomials b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff"); a = 1 - poly([x(nb+1:nb+na)]',"q","coeff") aSize = coeff(a);bSize = coeff(b) @@ -1,14 +1,50 @@ -// Estimates Discrete time BJ 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_bj,opt_err,resid] = bj(varargin) function sys = bj(varargin) +// Parameters Estimation of BJ(Box-Jenkins) model using Input Output time-domain data +// +// Calling Sequence +// sys = bj(ioData,[nb nc nd nf nk]) +// +// Parameters +// ioData : iddata or [outputData inputData] ,matrix of nx2 dimensions, type plant data +// nb : non-negative integer number specified as order of the polynomial B(z^-1)+1 +// nc : non-negative integer number specified as order of the polynomial C(z^-1) +// nd : non-negative integer number specified as order of the polynomial D(z^-1) +// nf : non-negative integer number specified as order of the polynomial f(z^-1) +// nk : non-negative integer number specified as input output delay, Default value is 1 +// sys : idpoly type polynomial have estimated coefficients of B(z^-1),C(z^-1),D(z^-1) and f(z^-1) polynomials +// +// Description +// Fit BJ model on given input output data +// The mathematical equation of the BJ model +// <latex> +// begin{eqnarray} +// y(n) = \frac {B(q)}{D(q)}u(n) + \frac {C(q)}{D(q)}e(t) +// end{eqnarray} +// </latex> +// It is SISO type model. It minimizes the sum of the squares of nonlinear functions using Levenberg-Marquardt algorithm. +// sys ,an idpoly type class, have different fields that contains estimated coefficients, sampling time, time unit and other estimated data in Report object. +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.5];b = [0 2 3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// ioData = iddata(y,u,0.1) +// sys = bj(ioData,[2,2,2,2,1]) +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.5];b = [0 2 3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// ioData = [y,u] +// sys = bj(ioData,[2,2,2,2,1]) +// +// Authors +// Ashutosh Kumar Bhargava, Harpreet,Inderpreet + [lhs , rhs] = argn(); if ( rhs < 2 ) then errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "bj", rhs); @@ -51,12 +87,12 @@ function sys = bj(varargin) nk = n(5); end - // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + // 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 - _objfun(UDATA,p,nd,nc,nf,nb,nk); + e = YDATA - _objfunbj(UDATA,p,nd,nc,nf,nb,nk); endfunction tempSum = nb+nc+nd+nf p0 = linspace(0.5,0.9,tempSum)'; @@ -71,12 +107,12 @@ function sys = bj(varargin) f = poly([1; var(nb+nd+nc+1:nd+nc+nf+nb)]',"q","coeff"); t = idpoly(1,coeff(b),coeff(c),coeff(d),coeff(f),Ts) - // estimating the other parameters + // estimating the other parameters [temp1,temp2,temp3] = predict(z,t) [temp11,temp22,temp33] = pe(z,t) estData = calModelPara(temp1,temp11,sum(n(1:4))) - //pause + // pause t.Report.Fit.MSE = estData.MSE t.Report.Fit.FPE = estData.FPE t.Report.Fit.FitPer = estData.FitPer @@ -88,11 +124,11 @@ function sys = bj(varargin) sys = t endfunction -function yhat = _objfun(UDATA,x,nd,nc,nf,nb,nk) +function yhat = _objfunbj(UDATA,x,nd,nc,nf,nb,nk) x=x(:) q = poly(0,'q') tempSum = nb+nc+nd+nf - // making polynomials + // making polynomials b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff"); c = poly([1; x(nb+1:nb+nc)]',"q","coeff"); d = poly([1; x(nb+nc+1:nb+nc+nd)]',"q","coeff"); diff --git a/buildmacros.sce b/buildmacros.sce new file mode 100644 index 0000000..5a872bb --- /dev/null +++ b/buildmacros.sce @@ -0,0 +1,10 @@ +// This file is released under the 3-clause BSD license. See COPYING-BSD. + +function buildmacros() + macros_path = get_absolute_file_path("buildmacros.sce"); + tbx_build_macros(TOOLBOX_NAME, macros_path); +endfunction + +buildmacros(); +clear buildmacros; // remove buildmacros on stack + diff --git a/calModelPara.sci b/calModelPara.sci index e8fbc02..135cf06 100644 --- a/calModelPara.sci +++ b/calModelPara.sci @@ -11,6 +11,6 @@ function varargout = calModelPara(varargin) AICc = AIC * 2 * np * (np + 1)/(N - np - 1);//disp(AICc) nAIC = log(mse) + 2 * np/N;//disp(nAIC) BIC = N * log(mse) + N * size(y,'c') * (log(2 * %pi) + 1) + np * log(N);//disp(BIC) - //pause + // pause varargout(1) = struct('MSE',mse,'FPE',fpe,'FitPer',nrmse*100,'AIC',AIC,'AICc',AICc,'nAIC',nAIC,'BIC',BIC) endfunction diff --git a/cleanmacros.sce b/cleanmacros.sce new file mode 100644 index 0000000..a6428b2 --- /dev/null +++ b/cleanmacros.sce @@ -0,0 +1,20 @@ +// ==================================================================== +// This file is released under the 3-clause BSD license. See COPYING-BSD. +// ==================================================================== +function cleanmacros() + + libpath = get_absolute_file_path("cleanmacros.sce"); + + binfiles = ls(libpath+"/*.bin"); + for i = 1:size(binfiles,"*") + mdelete(binfiles(i)); + end + + mdelete(libpath+"/names"); + mdelete(libpath+"/lib"); +endfunction + +cleanmacros(); +clear cleanmacros; // remove cleanmacros on stack + +// ==================================================================== diff --git a/compare.sci b/compare.sci index 4723130..f171863 100644 --- a/compare.sci +++ b/compare.sci @@ -1,145 +1,143 @@ + function varargout = compare(varargin) + +// Make comparision between plant output and estimated output +// +// Calling Sequence +// compare(plantData,sys) +// [yData,tData,fData] = compare(plantData,sys) +// Parameters +// plantData : iddata type or nx2 matrix +// sys : idpoly type polynomial +// yData : one step ahead estimated output response +// tData : time series data +// fData : real number represent the normalized root mean square (NRMS) of the estimated output response +// Description +// compare function do the estimation of the one step ahead output response of the iddata type polynomial and compare it with the actual output response data with (NRMS). +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// y = sim(u,sys)+rand(1024,1) +// plantData = iddata(y,u,0.1) +// compare(plantData,sys) +// Authors +// Ashutosh Kumar Bhargava + [lhs,rhs] = argn(0) -//------------------------------------------------------------------------------ -// checking the number of inputs - if rhs < 2 || rhs > 3 then - error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"compare")) + sampleData = varargin(1) + sysData = varargin(2) + yData = [] + uData = [] + Ts = 0 + if lhs > 3 then + error(msprintf(gettext("%s:Wrong number of output arguments.\n"),"compare")) end -//------------------------------------------------------------------------------ - data = varargin(1) - model = varargin(2) - if rhs == 3 then - kStep = varargin(3) - elseif rhs == 2 then - kStep = %inf + if type(varargin(1)) == 0 then + error(msprintf(gettext("%s:Null sample data matrix.\n"),"compare")) + end + // if the plant data stores in "iddata" function + if typeof(sampleData) == 'iddata' then + yData = sampleData.OutputData + uData = sampleData.InputData + Ts = sampleData.Ts + + // if the plant data stores in a matrix of nx2 -->([yData,uData]) + elseif typeof(sampleData) == 'constant' then + if size(sampleData,'c') ~= 2 then + error(msprintf(gettext("%s:Incorrect number of columns in sample data.\n"),"compare")) + end + yData = sampleData(:,1) + uData = sampleData(:,2) + Ts = 1 + else + error(msprintf(gettext("%s:Improper sample datas.\n"),"compare")) end -//------------------------------------------------------------------------------ - -// k step analysis - if typeof(kStep) <> 'constant' || isnan(kStep) then - error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"compare")) + if ~size(yData,'*') then + error(msprintf(gettext("%s:Sample data must contain plant output datas.\n"),"compare")) end -// if given k step is infinity or [] -// if isinf(kStep) || ~size(kStep,'*') then -// kStep = 1 -// end -// checking the dimensions -// if size(kStep,'*') <> 1 || (ceil(kStep)-kStep) then -// error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"predict")) -// end -//------------------------------------------------------------------------------ -// checking the plant model - if typeof(model) ~= 'idpoly' then - error(msprintf(gettext("%s:Plant model must be ""idpoly"" type.\n"),"compare")) + + if ~size(uData,'*') then + error(msprintf(gettext("%s:Sample data must contain plant input datas.\n"),"compare")) end - modelSampleTime = model.Ts - modelTimeUnit = model.TimeUnit -//------------------------------------------------------------------------------ -//checking the data type - if typeof(data) <> 'iddata' && typeof(data) <> 'constant' then - error(msprintf(gettext("%s:Sample data must be ""iddata"" type or ""n x 2"" matrix type.\n"),"compare")) + + if typeof(sysData) == 'state-space' then + elseif typeof(sysData) == 'rational' then + sysData = tf2ss(sysData) + elseif typeof(sysData) == 'idpoly' then + else + error(msprintf(gettext("%s: Wrong type for input argument \n#%d: State-space \n#%d: Transfer function \n#%d: System identification \n expected.\n"),"compare",1,2,3)) end -// checking the plant data - if typeof(data) == 'iddata' then - if ~size(data.OutputData,'*') || ~size(data.InputData,'*') then - error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"compare")) - end - plantSampleTime = data.Ts - plantTimeUnit = data.TimeUnit - data = [data.OutputData data.InputData] - //disp('iddata') - elseif typeof(data) == 'constant' then - if size(data,'c') ~= 2 then - error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"compare")) + sampleLength = size(yData,'r') + x0 = []// initial state + // for state space type systems + if typeof(sysData) == 'state-space' then +// sampleLength = size(yData,'r') + sysData = dscr(sysData,Ts) + ySys = [0] + for ii = 2:sampleLength + tempData = 0 + for jj = 1:ii-1 + tempData = tempData + (sysData.c)*(sysData.a)^(ii-jj-1)*(sysData.b)*uData(jj) + end + ySys = [ySys; tempData] end - plantSampleTime = model.Ts - plantTimeUnit = model.TimeUnit + x0 = zeros(size(sysData.A,'r'),1) + + // for system identification type system(OE,BJ system Models) + elseif typeof(sysData) == 'idpoly' then + Ts = sysData.Ts + zd = [yData uData] + ySys = compareBJ(sysData,zd) + // oe model comparision +// if size(sysData.a,'*') == 1 & size(sysData.c,'*') == 1 & size(sysData.d,'*') == 1 & size(sysData.b,'*') ~= 1 & size(sysData.f,'*') ~= 1 then +// ySys = compareOE(sysData,zd) +// else +// ySys = compareBJ(sysData,zd) +// end end -//------------------------------------------------------------------------------ -// comparing the sampling time - if modelSampleTime-plantSampleTime <> 0 then - error(msprintf(gettext("%s:The sample time of the model and plant data must be equal.\n"),"compare")) + if isrow(yData) then + yData = yData' end -// Comparing the time units - if ~strcmp(modelTimeUnit,plantTimeUnit) then - else - error(msprintf(gettext("%s:Time unit of the model and plant data must be equal.\n"),"compare")) + if isrow(ySys) then + ySys = ySys' end -//------------------------------------------------------------------------------ -// ckecking the k step size. if it greater than number of sample size then the -// k step will become 1 -// if kStep >= size(data,'r') then -// kStep = 1 -// end -//------------------------------------------------------------------------------ - //storing the plant data - // B(z) C(z) - // y(n) = ---------- u(n) + ---------- e(n) - // A(z)*F(z) A(z)*D(z) - aPoly = poly(model.a,'q','coeff'); - bPoly = poly(model.b,'q','coeff'); - cPoly = poly(model.c,'q','coeff'); - dPoly = poly(model.d,'q','coeff'); - fPoly = poly(model.f,'q','coeff'); - Gq = bPoly/(aPoly*fPoly) - Hq = cPoly/(aPoly*dPoly) - //disp(kStep) - if kStep == %inf then - //disp('in inf') - outData = sim(data(:,2),model) - else - if kStep == 1 then - Wkq = Hq^-1 - elseif kStep > 1 then - adCoeff = coeff(aPoly*dPoly);adCoeff = adCoeff(length(adCoeff):-1:1); - adPoly = poly(adCoeff,'q','coeff') - cCoeff = model.c;cCoeff = cCoeff(length(cCoeff):-1:1); - cPoly = poly(cCoeff,'q','coeff') - hBar = clean((ldiv(cPoly,adPoly,kStep))',0.00001) - hBarPoly = poly(hBar,'q','coeff') - Wkq = hBarPoly*Hq^-1//*bPoly/(dPoly*fPoly) - end - // pause - WkqGq = Wkq * Gq - tempWkqGq = coeff(WkqGq.den) - if tempWkqGq(1) <> 1 then - WkqGq.num = WkqGq.num/tempWkqGq(1) - WkqGq.den = WkqGq.den/tempWkqGq(1) - end - Wkq1 = 1-Wkq - tempWkq1 = coeff(Wkq1.den) - if tempWkq1(1) == 1 then - Wkq1.num = Wkq1.num/tempWkq1(1) - Wkq1.den = Wkq1.den/tempWkq1(1) + fit = fitValue(yData,ySys) + + if lhs == 1 then + varargout(1) = fit// zeros(size(sysData.A,'r'),1) + dataTime = 0:Ts:(sampleLength-1)*Ts + plot(dataTime',ySys) + plot(dataTime',yData,'m') + // plot(ySys) + // plot(yData,'m') + xData = 'Time(' + if typeof(sampleData) == 'constant' then + xData = xData + 'seconds)' + elseif typeof(sampleData) == 'iddata' then + xData = xData + sampleData.TimeUnit+')' end - z = poly(0,'z'); - WkqGqNum = horner(WkqGq.num,1/z);WkqGqDen = horner(WkqGq.den,1/z); - uPoly = WkqGqNum/WkqGqDen; - uPoly = syslin('d',uPoly);uData = flts(data(:,2)',uPoly); - - Wkq1Num = horner(Wkq1.num,1/z);Wkq1Den = horner(Wkq1.den,1/z); - yPoly = Wkq1Num/Wkq1Den; - yPoly = syslin('d',yPoly);yData = flts(data(:,1)',yPoly); - outData = (uData+yData)' - end - - tData = (1:size(data,'r'))'*plantSampleTime - fitData = fitValue(data(:,1),outData) - if lhs == 1 then - varargout(1) = [] - plot(tData,data(:,1),'m') - plot(tData,outData,'b') - legend('Plant Data','Model Data : '+string(fitData)+'%') - xtitle('Comparison of Time Response','Time ('+ plantTimeUnit+')','Amplitude') + h = gcf() + h.figure_name = 'Simulated Compare Plot' xgrid() + xtitle('Comparison of Time Response',xData,'Amplitude') + plant = 'plant : ' + string(fit)+ '% fit' + legend(['sys',plant]) elseif lhs == 2 then - varargout(1) = fitData - varargout(2) = outData + varargout(2) = fit + varargout(1) = x0 elseif lhs == 3 then - varargout(1) = outData - varargout(2) = tData - varargout(3) = fitData + varargout(3) = fit// zeros(size(sysData.A,'r'),1) + varargout(2) = x0 + yOutput = iddata(ySys,[],Ts) + if typeof(varargin(1)) == 'constant' then + yOutput.TimeUnit = 'seconds' + elseif typeof(varargin(1)) == 'iddata' then + yOutout.TimeUnit = sampleData.TimeUnit + end + varargout(1) = yOutput end endfunction + diff --git a/dataSlice.sci b/dataSlice.sci index 6b3c160..408bce9 100644 --- a/dataSlice.sci +++ b/dataSlice.sci @@ -1,26 +1,49 @@ + function sys = dataSlice(data,Start,End,Freq) +// Select sample data from iddata +// +// Calling Sequence +// h = dataSlice(plantData,Start,End,Ts) +// Parameters +// data : iddata type +// Start : non-negative integer index +// End : non-negative integer index, always greater than Start index +// Ts : sampling frequency, default value is 1 +// Description +// Extracts the samples in between Start and End index of the plant time series data,iddata type. For specified sampling frequency, it resamples the extracted data. +// +// Examples +// a = [1 0.5];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// y = sim(u,sys)+rand(1024,1) +// plantData = iddata(y,u,0.1) +// h = dataSlice(plantData,1,20,0.1) +// Authors +// Ashutosh Kumar Bhargava,Bhushan Manjarekar + [lhs,rhs] = argn() - // storing the model data + // storing the model data modelData = data - // storing the statrting point + // storing the statrting point try startData = Start catch startData = 1 end - //storing the end point + // storing the end point try endData = End catch endData = LastIndex(data) end - //Storing the frequency + // Storing the frequency try freqData = Freq catch freqData = 1 end - // error message generate + // error message generate if startData > endData then error(msprintf(gettext("%s:Start index can not greater than End index.\n"),"dataSlice")) end @@ -33,13 +56,13 @@ function sys = dataSlice(data,Start,End,Freq) if ~freqData || size(freqData,'*') ~= 1 then error(msprintf(gettext("%s:Frequency must be non negative scalar number.\n"),"dataSlice")) end - //-------------------------------------------------------------------------- + // -------------------------------------------------------------------------- if typeof(modelData) == 'constant' then Ts = 1 elseif typeof(modelData) == 'iddata' then Ts = modelData.Ts end - //-------------------------------------------------------------------------- + // -------------------------------------------------------------------------- if freqData> Ts || modulo(Ts,freqData) then warning(msprintf(gettext("%s:inconsistent frequency.\n"),"dataSlice")) freqData = Ts @@ -48,7 +71,16 @@ function sys = dataSlice(data,Start,End,Freq) temp = modelData(startData:Ts/freqData:endData,:) elseif typeof(modelData) == 'iddata' then tempY = modelData.OutputData;tempU = modelData.InputData - tempY = tempY(startData:Ts/freqData:endData,:);tempU = tempU(startData:Ts/freqData:endData,:) + if ~size(tempY,'r') then + tempY = [] + else + tempY = tempY(startData:Ts/freqData:endData,:); + end + if ~size(tempU,'r') then + tempU = [] + else + tempU = tempU(startData:Ts/freqData:endData,:) + end temp = iddata(tempY,tempU,Ts/freqData) temp.TimeUnit = modelData.TimeUnit end @@ -56,12 +88,12 @@ function sys = dataSlice(data,Start,End,Freq) endfunction function varargout = LastIndex(modelData) - //finding the sample size + // finding the sample size if typeof(modelData) == "constant" then varargout(1) = length(modelData(:,1)) elseif typeof(modelData) == "iddata" then - temp = modelData.OutputData - varargout(1) = length(temp(:,1)) + temp = max(size(modelData.OutputData,'r'),size(modelData.InputData,'r')) + varargout(1) = temp end endfunction diff --git a/estpoly.sci b/estpoly.sci index f1ea48d..4fe9362 100644 --- a/estpoly.sci +++ b/estpoly.sci @@ -1,14 +1,53 @@ -// 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) - + +// Parameters Estimation of generalized discrete time model using Input Output time-domain data +// +// Calling Sequence +// sys = bj(ioData,[na nb nc nd nf nk]) +// +// Parameters +// ioData : iddata or [outputData inputData] ,matrix of nx2 dimensions, type plant data +// na : non-negative integer number specified as order of the polynomial A(z^-1) +// nb : non-negative integer number specified as order of the polynomial B(z^-1)+1 +// nc : non-negative integer number specified as order of the polynomial C(z^-1) +// nd : non-negative integer number specified as order of the polynomial D(z^-1) +// nf : non-negative integer number specified as order of the polynomial F(z^-1) +// nk : non-negative integer number specified as input output delay, Default value is 1 +// sys : idpoly type polynomial have estimated coefficients of A(z^-1), B(z^-1),C(z^-1),D(z^-1) and F(z^-1) polynomials +// +// Description +// Fit generalized discrete time model on given input output data +// The mathematical equation of the generalized discrete time model +// <latex> +// begin{eqnarray} +// A(q)y(n) = \frac {B(q)}{D(q)}u(n) + \frac {C(q)}{D(q)}e(t) +// end{eqnarray} +// </latex> +// It is SISO type model. It minimizes the sum of the squares of nonlinear functions using Levenberg-Marquardt algorithm. +// +// sys ,an idpoly type class, have different fields that contains estimated coefficients, sampling time, time unit and other estimated data in Report object. +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.2];b = [0 0.2 0.3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// ioData = iddata(y,u,0.1) +// sys = estpoly(ioData,[2,2,2,2,2,1]) +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.2];b = [0 0.2 0.3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// ioData = [y,u] +// sys = estpoly(ioData,[2,2,2,2,2,1]) +// +// Authors +// Ashutosh Kumar Bhargava + [lhs , rhs] = argn(); if ( rhs < 2 ) then errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "estpoly", rhs); @@ -51,26 +90,26 @@ function sys = estpoly(varargin) nk = n(6); end - // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + // 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); + 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) + // disp(errl) err = (norm(errl)^2); - //disp(err) + // 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"); +// 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"); @@ -78,14 +117,14 @@ function sys = estpoly(varargin) 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) + // t = sys;// idpoly(1,coeff(b),coeff(c),coeff(d),coeff(f),Ts) - // estimating the other parameters + // estimating the other parameters [temp1,temp2,temp3] = predict(z,t) [temp11,temp22,temp33] = pe(z,t) - //pause + // pause estData = calModelPara(temp1,temp11,na+nb+nc+nd+nf) - //pause + // pause t.Report.Fit.MSE = estData.MSE t.Report.Fit.FPE = estData.FPE t.Report.Fit.FitPer = estData.FitPer @@ -96,14 +135,14 @@ function sys = estpoly(varargin) t.TimeUnit = unit sys = t - //sys.TimeUnit = unit + // sys.TimeUnit = unit endfunction -function yhat = _oestpolyfun(UDATA,x,na,nb,nc,nd,nf,nk)//(UDATA,x,nd,nc,nf,nb,nk) +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 + // 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"); @@ -113,7 +152,7 @@ function yhat = _oestpolyfun(UDATA,x,na,nb,nc,nd,nf,nk)//(UDATA,x,nd,nc,nf,nb,nk if size(bd,"*") == 1 then bd = repmat(0,nb+nd+1,1) end - //pause + // pause maxDelay = max([length(bd) length(cf) length(fc_d)]) yhat = [YDATA(1:maxDelay)] for k=maxDelay+1:size(UDATA,"*") @@ -130,7 +169,7 @@ function yhat = _oestpolyfun(UDATA,x,na,nb,nc,nd,nf,nk)//(UDATA,x,nd,nc,nf,nb,nk for i = 2:size(cf,"*") cfadd = cfadd + cf(i)*yhat(k-i+1) end - //pause + // pause yhat = [yhat; [ bdadd + fc_dadd - cfadd ]]; end endfunction @@ -1,18 +1,50 @@ function varargout = etfe(varargin) + +// Estimate empirical transfer function +// +// Calling Sequence +// frdData = etfe(plantData,n) +// Parameters +// plantData : iddata type +// n : frequency sample spacing, default value is 128 +// frdData : frd type object +// Description +// etfe function takes time domain plant data,iddata type, and estimate the empirical transfer function by taking the ratio of the fourier transforms +// of the output and the input variables +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// y = sim(u,sys)+rand(1024,1) +// plantData = iddata(y,u,0.1) +// frdData = etfe(plantData) +// Authors +// Ashutosh Kumar Bhargava,Bhushan Manjarekar + [lhs,rhs] = argn() + if rhs > 2 || rhs < 1 then + error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"etfe")) + end data = varargin(1) + if typeof(data) <> "iddata" then + error(msprintf(gettext("%s:Plant time series data must be ""iddata"" type.\n"),"etfe")) + end if rhs == 1 then n = 128 elseif rhs == 2 then n = varargin(2) end y = data.OutputData; - u = data.InputData + u = data.InputData; + if ~size(u,'r') then + error(msprintf(gettext("%s:Non zero input data point needed.\n"),"etfe")) + elseif ~size(y,'r') then + error(msprintf(gettext("%s:Non zero output data point needed.\n"),"etfe")) + end N = size(y,'r') - v = linspace(1,N,n) - y1 = y(v);u1 = u(v) - //y1 = y((1:(N-1)/(n-1):N));u1 = u((1:(N-1)/(n-1):N)) - //y1($) = y(N);u1($) = u(N) + + y1 = y((1:ceil((N-1)/(n-1)):N));u1 = u((1:ceil((N-1)/(n-1)):N)) + y1($) = y(N);u1($) = u(N) data12 = [y1,u1] z = [fft(y1),fft(u1)] z = z/size(z,'r') @@ -24,6 +56,6 @@ function varargout = etfe(varargin) output = data(:,1)+%i*data(:,2) resp = output(1:ceil(length(output)/2)) frq = (1: ceil(n/2)) * %pi/floor(n/2) - output = frd(frq,resp,1) + output = frd(frq',resp,1) varargout(1)= output endfunction diff --git a/fitValue.sci b/fitValue.sci new file mode 100644 index 0000000..75db2a6 --- /dev/null +++ b/fitValue.sci @@ -0,0 +1,14 @@ +function varargout = fitValue(varargin) + yData = varargin(1) + ySys = varargin(2) + fitValueNum = 0 + fitValueDen = 0 + sampleLength = size(yData,'r') + yDataMean = mean(yData) + for ii = 1:sampleLength + fitValueNum = fitValueNum + (yData(ii)-ySys(ii))^2 + fitValueDen = fitValueDen + (yData(ii)-yDataMean)^2 + end + fit = 100*(1-sqrt(fitValueNum/fitValueDen)) + varargout(1) = fit +endfunction @@ -1,4 +1,32 @@ function varargout = fitch(varargin) + +// Characteristics of estimated model +// +// Calling Sequence +// estData = fitch(sys) +// Parameters +// sys : idpoly type polynomial +// estData : stuct type variable have following objects +// MSE : Mean Square Error, show the wellness of the estimated model on plant data +// FPE : Final Prediction Error +// FitPer : Normalized root mean squared error (NRMSE), show the precentage fit of the estimated model on plant data +// AIC : Raw Akaike Information Citeria (AIC), show the model quality +// AICc : Small sample-size corrected AIC +// nAIC : Normalized AIC +// BIC : Bayesian Information Criteria (BIC) +// Description +// fitch function represent other calculated model characteristics to show the wellness of the model on different scales. +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys1 = idpoly(a,b,'Ts',0.1) +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// y = sim(u,sys)+rand(1024,1) +// plantData = iddata(y,u,0.1) +// sys = arx(plantData,[2,2,1]) +// estData = fitch(sys) +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + [lhs , rhs] = argn() if ( rhs <> 1 ) then errmsg = msprintf(gettext("%s: Wrong number of input arguments"), "fitch"); @@ -6,15 +34,15 @@ function varargout = fitch(varargin) elseif typeof(varargin(1)) <> "idpoly" then error(msprintf(gettext("%s:Input model must be ""idpoly"" type.\n"),"fitch")) end - model = varargin(1) - MSE = model.Report.Fit.MSE - FPE = model.Report.Fit.FPE + model = varargin(1) + MSE = model.Report.Fit.MSE + FPE = model.Report.Fit.FPE FitPer = model.Report.Fit.FitPer - AIC = model.Report.Fit.AIC - AICc = model.Report.Fit.AICc - nAIC = model.Report.Fit.nAIC - BIC = model.Report.Fit.BIC - t = tlist(['fitch','MSE','FPE','FitPer','AIC','AICc','nAIC','BIC'],MSE,FPE,FitPer,AIC,AICc,nAIC,BIC) + AIC = model.Report.Fit.AIC + AICc = model.Report.Fit.AICc + nAIC = model.Report.Fit.nAIC + BIC = model.Report.Fit.BIC + t = tlist(['fitch','MSE','FPE','FitPer','AIC','AICc','nAIC','BIC'],MSE,FPE,FitPer,AIC,AICc,nAIC,BIC) varargout(1) = t endfunction @@ -1,17 +1,40 @@ function varargout = frd(varargin) + + +// Stores frequency and response data +// +// Calling Sequence +// plantData = frd(respData,frdData,Ts) +// Parameters +// frdData : nx1 matrix of non-decreasing frequency points +// respData : nx1 matrix of the frequency response +// Ts : non-negative real number +// plantData : frd type module +// Description +// It is a frd type module that stores the frequency and response data with sampling time Ts. The time unite is in second and frequency unit rad/sec. +// Examples +// frdData = 0:1024 +// frdData=frdData'; +// respData = rand(1024,1) +// Ts = 0.1 +// plantData = frd(frdData,respData,Ts) +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + + [lhs,rhs] = argn(0) if rhs < 2 || rhs > 4 then errmsg = msprintf(gettext("%s: Wrong numbers of input arguments."), "frd"); error(errmsg) end - frequency = varargin(2) + frequency = varargin(1) freqUnit = 'rad/TimeUnit' if ~iscolumn(frequency) then - errmsg = msprintf(gettext("%s: frequency must be a finite row vector."), "frd"); + errmsg = msprintf(gettext("%s: frequency must be a finite column vector."), "frd"); error(errmsg) end - respData = varargin(1) -// pause + respData = varargin(2) +// pause if size(frequency,'r') <> size(respData,'r') then errmsg = msprintf(gettext("%s: input output matrix dimension must be equal."), "frd"); error(errmsg) @@ -25,13 +48,13 @@ function varargout = frd(varargin) errmsg = msprintf(gettext("%s: Sampling time must be a scalar non negative real number."), "frd"); error(errmsg) end - // saving the spectrum value + // saving the spectrum value if rhs == 4 then spect = varargin(4) else spect = [] end - /// matching its dimensions + // / matching its dimensions if ~size(spect) then elseif size(frequency,'r') <> size(spect,'r') then errmsg = msprintf(gettext("%s: Numbers of power spectra must be equal to the numbers of frequency."), "frd"); @@ -42,7 +65,7 @@ function varargout = frd(varargin) varargout(1) = t endfunction -//overloading +// overloading function %frd_p(varargin) myTlist = varargin(1) f = fieldnames(myTlist) @@ -62,7 +85,7 @@ function %frd_p(varargin) temp = temp + ' ' end temp = temp + string(real(respData(ii))) - // temp = string(real(respData(ii))) + // temp = string(real(respData(ii))) if imag(respData(ii)) > 0 then temp = temp +"+" end @@ -70,8 +93,8 @@ function %frd_p(varargin) else temp = temp + string(imag(respData(ii))) +"i" end - // temp = temp + string(imag(respData(ii))) + " i" - mprintf("\n\t %f \t %s",freqData(ii),temp)//real(respData(ii)),imag(respData(ii))) + // temp = temp + string(imag(respData(ii))) + " i" + mprintf("\n\t %f \t %s",freqData(ii),temp)// real(respData(ii)),imag(respData(ii))) end mprintf("\n\n") end diff --git a/frdplot.sci b/frdplot.sci index c923429..f091e54 100644 --- a/frdplot.sci +++ b/frdplot.sci @@ -1,4 +1,23 @@ function frdplot(varargin) + +// Plot the frequency domain response +// +// Calling Sequence +// frdplot(plantData) +// Parameters +// frdData : frd type module +// +// Description +// plantData must be frd type. Function takes the frequency and response data, and plot the bode plot. +// Examples +// frdData = (1:100)'; +// respData = rand(100,1) + %i * rand(100,1); +// plantData = frd(frdData,respData); +// frdplot(plantData) +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + + [lhs rhs] = argn(0) if rhs <> 1 then error(msprintf(gettext("%s: Wrong number of input arguments."),"frdplot")) @@ -1,4 +1,35 @@ function varargout = iddata(varargin)
+
+// Stores plant input-output data
+//
+// Calling Sequence
+// plantData = iddata(yData,uData,Ts)
+// plantData = iddata([],uData,Ts)
+// plantData = iddata(yData,[],Ts)
+// plantData = iddata(yData,uData)
+// plantData = iddata(yData)
+//
+// Parameters
+// uData : nx1 matrix of the plant input data
+// yData : nx1 matrix of the plant output data
+// Ts : non-negative real number
+// plantData : iddata type module
+//
+// Description
+// It is a iddata type module that stores the plant input and output data with sampling time Ts. The time unite is in second.
+//
+// Examples
+// uData = idinput(1024,'PRBS',[0 1/20],[-1 1])
+// yData = rand(1024,1)
+// Ts = 0.1
+// plantData1 = iddata(yData,uData,Ts)
+// plantData2 = iddata(yData,[],Ts)
+// plantData3 = iddata([],uData,Ts)
+//
+// Authors
+// Ashutosh Kumar Bhargava
+
+
[lhs,rhs] = argn(0)
if rhs == 2 | rhs == 1 then
Ts = 1
diff --git a/iddataplot.sci b/iddataplot.sci index 4b1818d..8261c9c 100644 --- a/iddataplot.sci +++ b/iddataplot.sci @@ -1,4 +1,30 @@ function iddataplot(varargin) + +// Plot the iddata class +// +// Calling Sequence +// iddataplot(plantData) +// Parameters +// plantData : iddata type module +// +// Description +// Takes the iddata type input. Plot the respective input and output data +// +// Examples +// uData = idinput(1024,'PRBS',[0 1/20],[-1 1]); +// yData = rand(1024,1); +// Ts = 0.1; +// plantData1 = iddata(yData,uData,Ts) +// iddataplot(plantData1) +// plantData2 = iddata(yData,[],Ts) +// iddataplot(plantData2) +// plantData3 = iddata([],uData,Ts) +// iddataplot(plantData3) +// +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + + [lhs rhs] = argn(0) if rhs <> 1 then error(msprintf(gettext("%s: Wrong number of input arguments."),"iddataplot")) @@ -7,19 +33,19 @@ function iddataplot(varargin) if typeof(iddataData) <> 'iddata' then error(msprintf(gettext("%s:Wrong type for input argument %d: ""iddata"" expected."),"iddataplot",1)) end -// figure() -// xtitle('Plant i/o Data','Time('+iddataData.TimeUnit+')','Amplitude') +// figure() +// xtitle('Plant i/o Data','Time('+iddataData.TimeUnit+')','Amplitude') uData = iddataData.InputData yData = iddataData.OutputData timeLength = max(size(uData,'r'),size(yData,'r')) timeData = ((1:1:timeLength)*iddataData.Ts)' - //timeData = timeData(1:length(iddataData)-1) + // timeData = timeData(1:length(iddataData)-1) if size(uData,'*') && size(yData,'*') then firstIndex = 2 else firstIndex = 1 end - // ploting y data + // ploting y data h = gcf() if h.figure_name == "Plant Input Output Data"; then clf() @@ -34,7 +60,7 @@ function iddataplot(varargin) outputString = [outputString 'y'+string(ii)] end end - //disp(outputString) + // disp(outputString) for ii = 1:secondIndex subplot(firstIndex,secondIndex,ii);plot(timeData,yData(:,ii)) xtitle(outputString(ii)) @@ -42,7 +68,7 @@ function iddataplot(varargin) end if size(uData,'*') then secondIndex = size(uData,'c') - // disp(secondIndex) + // disp(secondIndex) if secondIndex == 1 then outputString = 'u' else @@ -51,7 +77,7 @@ function iddataplot(varargin) outputString = [outputString 'u'+string(ii)] end end - //disp(outputString) + // disp(outputString) for ii = 1:secondIndex if size(yData,'*') then temp = 1 @@ -59,11 +85,11 @@ function iddataplot(varargin) temp = 0 end subplot(firstIndex,secondIndex,ii+secondIndex*temp);plot2d(timeData,uData(:,ii),2) - //pause + // pause xtitle(outputString(ii)) end end h = gcf() - //disp(h) + // disp(h) h.figure_name= "Plant Input Output Data"; endfunction diff --git a/identTime.sci b/identTime.sci index c4c47c8..1395899 100644 --- a/identTime.sci +++ b/identTime.sci @@ -1,4 +1,25 @@ function varargout = identTime(varargin) + +// Returns the time series data +// +// Calling Sequence +// timeData = identTime(plantData) +// Parameters +// plantData : iddata type variable +// timeData : identTime type variable +// Description +// identTime function takes plantData,iddata type, and generate time vector data,frequency as number of samples per time unit,number of samples. The timeData is identType +// variable. +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys1 = idpoly(a,b,'Ts',0.1) +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// y = sim(u,sys)+rand(1024,1) +// plantData = iddata(y,u,0.1) +// timeData = identTime(plantData) +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + [lhs , rhs] = argn() if ( rhs <> 1 ) then errmsg = msprintf(gettext("%s: Wrong number of input arguments"), "identTime"); @@ -7,7 +28,7 @@ function varargout = identTime(varargin) error(msprintf(gettext("%s:Input model must be ""iddata"" type.\n"),"identTime")) end plantData = varargin(1) - //disp('yolo') + // disp('yolo') inputData = plantData.InputData;inputData = size(inputData,'r') outputData = plantData.OutputData;outputData = size(outputData,'r') sampleNumb = max(inputData,outputData) diff --git a/idinput.sci b/idinput.sci new file mode 100644 index 0000000..7fd62e6 --- /dev/null +++ b/idinput.sci @@ -0,0 +1,213 @@ +function[u]=idinput(N,types,band,levels) +// generates random binary input signal +// +// Calling Seqence +// u=idinput(N); +// u=idinput([n,nu]) +// u=idinput([n,nu,m]) +// u=idinput(__,type) +// u=idinput(__,type,band) +// u=idinput(__,type,band,levels) + +// +// Parameters +// N : no of binary signals +// [n,nu] : nu channel signal of length n. +// [n,nu,m]: nu channel signal of length nxm and period m. +// type : "rbs" ,"prbs" +// band : frequency band of the signal. +// levels : binary output values +// +// Description +// u = idinput(n) returns a single-channel random binary signal u of length N. By default it’s values are either -1 or 1. +// u = idinput([n,nu]) returns an nu-channel random binary input signal, where each channel signal has length n,the +// signals in each channel are independent from each other. +// u = idinput([n,nu,m]) returns an Nu-channel periodic random binary input signal with specified period and number of +// periods. Each input channel signal is of length m*n. +// u=idinput(___,type) specifies the type of input to be generated: + // 'rbs' — Random binary signal. + // 'prbs'—pseudorandom binary signal. +// u =idinput(___,type,band) specifies the frequency band of the signal. +// u =idinput(___,type,band,levels) specifies the level of the binary generated signal +// +// +// Examples +// u=idinput([20,2,2],'rbs',[0 0.3],[-1 2]); +// +// Author +// Ayush Kumar + + + [lhs,rhs]=argn(0) + if rhs<4 then // checking the input given with function + levels=[-1 1] // arguments and allocating default value if required + end + if rhs<3 then + band=[0 0.5] + end + if rhs<2 then + types='rbs' + end + if size(N,2)==3 then + P=N(1); + nu=N(2); + M=N(3); + elseif size(N,2)==2 then + P=N(1); + nu=N(2); + M=1; + elseif size(N,2)==1 then + P=N; + nu=1; + M=1; + else + error('erroininputargument'); + end + + if levels(2)<levels(1) then // error checking + error('levels (1)should be less then levels(2)'); + end + if type(types)~=10 then + error('types can be rbs only input valid argument'); + end + if(band(2)>1) then + error('band should be in b/w 0-1'); + end + + function[vec]=randnu(NN,nuu) // nuu channel random normally distributed signal generator + [lhs,rhs]=argn(0) + if rhs<2 then + nuu=1 + end + rand('seed',getdate('s')) + for j=1:1:nuu + for i= 1:NN + vec(i,j)=rand(NN,'normal') + end + end + endfunction + + + if convstr(types)=='rbs' then + u=randnu(5*P,nu); + band=band/2; + if ~and(band==[0 0.5]) then + if(band(1)==0) then + [hz]=iir(8,'lp','butt',[band(2) band(1)],[0 0]); // 8th order butterwoth filter + num=hz(2); + den=hz(3); + for i=1:1:nu + y(:,i)=filter(num,den,u(:,i)); + end + elseif(band(2)==0.5) then + [hz]=iir(8,'hp','butt',[band(1) band(2)],[0 0]); + num=hz(2); + den=hz(3); + for i=1:1:nu + y(:,i)=filter(num,den,u(:,i)); + end + else + [hz]=iir(8,'bp','butt',band,[0 0]); + num=hz(2); + den=hz(3); + for i=1:1:nu + y(:,i)=filter(num,den,u(:,i)); + end + end + + + u = sign(y(2*P+1:3*P,:)); // taking the middle terms + u = (levels(2)-levels(1))*(u+1)/2+levels(1);// adjusting binary values according to levels + else + u = sign(u(2*P+1:3*P,:)); + u = (levels(2)-levels(1))*(u+1)/2+levels(1); + + end + elseif convstr(types)=='prbs' then + clockP = floor(1/band(2)); + possP = 2.^(3:18)-1; + P1 = max(possP(P/clockP-possP>=0)); + if isempty(P1) + P1 = 7; + end + n = find(P1==possP)+2; + + if (clockP*P1~=P) + if M>1 + warning(msprintf(gettext("%s: period of prbs signal changed to %d and length change to %d"),"idinput",clockP*P1,P1*M*clockP)); + P = P1*clockP; + else + n = min(n+1,18); + P1 = 2^n -1; + warning(msprintf(gettext("%s :generated signal is first %d values of a sequence of length %d"),"idinput",P,clockP*P1)); + end + end + + P1 = 2^n-1; + + if n<3 || n>18, + error(msprintf(gettext('Bad conditioning'))); + end + fi = -ones(n,1); + if n==3 + ind = [1,3]; + elseif n==4 + ind = [1,4]; + elseif n==5 + ind = [2,5]; + elseif n==6 + ind = [1,6]; + elseif n==7 + ind = [1,7]; + elseif n==8 + ind = [1,2,7,8]; + elseif n==9 + ind = [4,9]; + elseif n==10 + ind = [3,10]; + elseif n==11 + ind = [9,11]; + elseif n==12 + ind = [6,8,11,12]; + elseif n==13 + ind = [9,10,12,13]; + elseif n==14 + ind = [4,8,13,14]; + elseif n==15 + ind = [14,15]; + elseif n==16 + ind = [4,13,15,16]; + elseif n==17 + ind = [14,17]; + elseif n==18 + ind = [11,18]; + end + for t = 1:clockP:P1*clockP + u(t:t+clockP-1,1) = ones(clockP,1)*fi(n); // LL%% multivariable !! + fi = [prod(fi(ind));fi(1:n-1,1)]; + end + u = (levels(2)-levels(1))*(u+1)/2+levels(1); + + u = u(1:P,1); + if nu >1 + u1 = [u;u]; + shift = floor(P/nu); + for ku = 2:nu + u = [u,u1(shift*(ku-1)+1:P+shift*(ku-1))]; + end + end + else + error(msprintf(gettext('type can be rbs or prbs only'))); + end + if M>1 then // generating periodic input if no. of periods>1 + uu = u; + for i = 2:M + u = [uu;u]; + end + end +endfunction + + + + + @@ -1,4 +1,34 @@ function sys = idpoly(varargin) + +// Stores the identification coefficients +// +// Calling Sequence +// sys = idpoly(aCoeff,bCoeff,cCoeff,dCoeff,fCoeff,Ts) +// +// Parameters +// aCoeff : 1xn coefficient matrix of the polynomial A(z^-1) +// bCoeff : 1xn coefficient matrix of the polynomial B(z^-1) +// cCoeff : 1xn coefficient matrix of the polynomial C(z^-1) +// dCoeff : 1xn coefficient matrix of the polynomial D(z^-1) +// fCoeff : 1xn coefficient matrix of the polynomial F(z^-1) +// +// Description +// It is a idpoly type polynomials that stores the identification coefficients A,B,C,D, and F with sampling time Ts. The time unite is in second. It stores other parameters +// in Report object. +// Examples +// aCoeff = [1 rand(1,2)] +// bCoeff = [0 rand(1,3)] +// cCoeff = [1 rand(1,2)] +// dCoeff = [1 rand(1,2)] +// fCoeff = [1 rand(1,2)] +// Ts = 0.1 +// sys1 = idpoly(aCoeff,bCoeff,cCoeff,dCoeff,fCoeff,Ts) +// sys2 = idpoly(aCoeff,bCoeff,cCoeff,'Ts',Ts) +// sys1 = idpoly(1,bCoeff,cCoeff,1,fCoeff,Ts) +// +// Authors +// Ashutosh Kumar Bhargava + [lhs,rhs] = argn(0) tempCell = cell(6,1) tempTs = 0 @@ -12,7 +42,7 @@ function sys = idpoly(varargin) tempCell{ii,1} = varargin(ii) end end - //disp(tempCell) + // disp(tempCell) for ii = 1:6 if ~size(cell2mat(tempCell(ii,1)),'*') & ii ~= 6 then tempCell{ii,1} = 1 @@ -20,10 +50,10 @@ function sys = idpoly(varargin) tempCell{ii,1} = -1 end end - //storing the data in A,B,C,D,F matrix - // B(z) D(z) - // y(n)=---------- u(n)+ ----------- e(n) - // A(z)F(z) A(z)D(z) + // storing the data in A,B,C,D,F matrix + // B(z) D(z) + // y(n)=---------- u(n)+ ----------- e(n) + // A(z)F(z) A(z)D(z) A = cell2mat(tempCell(1,1)); B = cell2mat(tempCell(2,1)); C = cell2mat(tempCell(3,1)); D = cell2mat(tempCell(4,1)); F = cell2mat(tempCell(5,1)); Ts = cell2mat(tempCell(6,1)); @@ -40,7 +70,7 @@ function sys = idpoly(varargin) errors = [0 0 0] report = struct('Fit',report,'Uncertainty',errors) t = tlist(['idpoly','a','b','c','d','f','Variable','TimeUnit','Ts','Report'],A,B,C,D,F,'z^-1','seconds',Ts,report) - //t = tlist(['idpoly','a','b','c','d','f','Variable','TimeUnit','Ts'],A,B,C,D,F,'z^-1','seconds',Ts) + // t = tlist(['idpoly','a','b','c','d','f','Variable','TimeUnit','Ts'],A,B,C,D,F,'z^-1','seconds',Ts) sys = t endfunction @@ -48,7 +78,7 @@ endfunction function %idpoly_p(mytlist) f = fieldnames(mytlist) - //A polynomial + // A polynomial if mytlist(f(1)) == 1 && size(mytlist(f(1)),'*') == 1 then else mprintf('\n A(z) =') @@ -56,7 +86,7 @@ function %idpoly_p(mytlist) mprintf('%s\n',temp) end - //B polynomial + // B polynomial if mytlist(f(2)) == 1 then else mprintf('\n B(z) =') @@ -64,14 +94,14 @@ function %idpoly_p(mytlist) mprintf('%s\n',temp) end - //C polynomial + // C polynomial if mytlist(f(3)) == 1 && size(mytlist(f(3)),'*') == 1 then else mprintf('\n C(z) =') temp = poly2str(mytlist(f(3))) mprintf('%s\n',temp) end - //D polynomial + // D polynomial if mytlist(f(4)) == 1 && size(mytlist(f(4)),'*') == 1 then elseif size(mytlist(f(4)),'*') > 1 then mprintf('\n D(z) =') @@ -79,7 +109,7 @@ function %idpoly_p(mytlist) mprintf('%s\n',temp) end - //F polynomial + // F polynomial if mytlist(f(5)) == 1 && size(mytlist(f(5)),'*') == 1 then else mprintf('\n F(z) =') @@ -98,7 +128,7 @@ function %idpoly_p(mytlist) mprintf('%0.4f %s',mytlist.Ts,mytlist.TimeUnit) end end - //disp(mytlist.Ts) + // disp(mytlist.Ts) mprintf('\n') if mytlist.Report.Fit.MSE then temp = ['MSE','FPE','FitPer','AIC','AICc','nAIC','BIC'] @@ -130,8 +160,8 @@ function strout = poly2str(h) temp = strsubst(temp,'*','') temp = strsubst(temp,'+',' + ') [ind which]=strindex(temp,'-') - //disp(ind) -// disp(which) + // disp(ind) +// disp(which) if ind(1,1) ~= 2 then temp = ' ' + temp elseif ind(1,1) == 2 then diff --git a/impulseest.sci b/impulseest.sci index 8817bcc..0522a18 100644 --- a/impulseest.sci +++ b/impulseest.sci @@ -1,6 +1,35 @@ function varargout = impulseest(varargin) + +// Estimate impulse response and plot of idpoly type model +// +// Calling Sequence +// impulseest(sys) +// impulseData = impulseest(sys,flag) +// +// Parameters +// sys : idpoly type polynomial +// flag : boolean type variable +// impulseData : stores impulse response if the flag value is true +// +// Description +// impulseest function estimate and plot the impulse response of idpoly type function. +// +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// impulseest(sys); +// +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// flag = %T +// impulseData = impulseest(sys,flag) +// +// Authors +// Ashutosh Kumar Bhargava + [lhs,rhs] = argn(0) - //checking the number of inputs + // checking the number of inputs if rhs > 2 then error(msprintf(gettext("%s: Unexpected number of input arguments "),"impulseest")) end @@ -8,7 +37,7 @@ function varargout = impulseest(varargin) if typeof(modelData) <> "idpoly" then error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"impulseest")) end - //adding noise + // adding noise if rhs == 2 then noiseFlag = varargin(2) if typeof(noiseFlag) <> 'boolean' then @@ -65,7 +94,7 @@ function varargout = impulseest(varargin) uData = [1 zeros(1,n)] yData = flts(uData,sys) timeData = (0:(n))*modelData.Ts -// pause +// pause if noiseFlag then varargout(1) = yData' else @@ -1,5 +1,50 @@ function varargout = iv(varargin) +// Parameters Estimation of IV model by arbitrary instrumental variable method +// +// Calling Sequence +// sys = iv(ioData,[na nb nk]) +// sys = iv(ioData,[na nb nk],instData) +// Parameters +// ioData : iddata or [outputData inputData] ,matrix of nx2 dimensions, type plant data +// na : non-negative integer number specified as order of the polynomial A(z^-1) +// nb : non-negative integer number specified as order of the polynomial B(z^-1)+1 +// nk : non-negative integer number specified as input output delay, Default value is 1 +// instData : arbitrary instrument variable matrix. The size of instriment variable must be equal to the size of outputData +// sys : idpoly type polynomial have estimated coefficients of A(z^-1) and B(z^-1) polynomials +// +// Description +// Fit IV model on given input output data +// The mathematical equation of the IV model +// <latex> +// begin{eqnarray} +// A(q)y(n) = B(q)u(n-k) + e(t) +// end{eqnarray} +// </latex> +// It is SISO type model. Instrument variable method is use to estimate the cofficients. If user does not provide the arbitrary instrument variable matrix +// then the program extracte it by using ARX method. +// +// sys ,an idpoly type class, have different fields that contains estimated coefficients, sampling time, time unit and other estimated data in Report object. +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.2];b = [0 0.2 0.3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// ioData = iddata(y,u,0.1) +// sys = arx(ioData,[2,2,1]) +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.2];b = [0 0.2 0.3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// ioData = [y,u] +// sys = iv(ioData,[2,2,1]) +// +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + [lhs , rhs] = argn(0); if ( rhs < 2 || rhs > 3) then errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2 or 3"), "iv", rhs); @@ -58,19 +103,19 @@ function varargout = iv(varargin) end phif = zeros(noOfSample,na+nb) psif = zeros(noOfSample,na+nb) - // arranging samples of y matrix + // arranging samples of y matrix for ii = 1:na phif(ii+1:ii+noOfSample,ii) = yData psif(ii+1:ii+noOfSample,ii) = x end - // arranging samples of u matrix + // arranging samples of u matrix for ii = 1:nb phif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData psif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData end lhs = psif'*phif lhsinv = pinv(lhs) - //pause + // pause theta = lhsinv * (psif)' * [yData;zeros(n,1)] ypred = (phif * theta) ypred = ypred(1:size(yData,'r')) @@ -79,12 +124,12 @@ function varargout = iv(varargin) vcov = sigma2 * pinv((phif)' * phif) t = idpoly([1; -theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,1) - // estimating the other parameters - [temp1,temp2,temp3] = predict(z,t) - [temp11,temp22,temp33] = pe(z,t) + // estimating the other parameters + [temp1,temp2,temp3] = predict(plantData,t) + [temp11,temp22,temp33] = pe(plantData,t) estData = calModelPara(temp1,temp11,na+nb) - //pause + // pause t.Report.Fit.MSE = estData.MSE t.Report.Fit.FPE = estData.FPE t.Report.Fit.FitPer = estData.FitPer @@ -93,6 +138,6 @@ function varargout = iv(varargin) t.Report.Fit.nAIC = estData.nAIC t.Report.Fit.BIC = estData.BIC t.TimeUnit = unit - //sys = t + // sys = t varargout(1) = t endfunction @@ -1,10 +1,53 @@ function varargout = iv4(varargin) +// Parameters Estimation of IV4 model by four stage instrumental variable method +// +// Calling Sequence +// sys = iv(ioData,[na nb nk]) +// +// Parameters +// ioData : iddata or [outputData inputData] ,matrix of nx2 dimensions, type plant data +// na : non-negative integer number specified as order of the polynomial A(z^-1) +// nb : non-negative integer number specified as order of the polynomial B(z^-1)+1 +// nk : non-negative integer number specified as input output delay, Default value is 1 +// sys : idpoly type polynomial have estimated coefficients of A(z^-1) and B(z^-1) polynomials +// +// Description +// Fit IV4 model on given input output data +// The structure of sys is ARX type.The mathematical equation is given here +// <latex> +// begin{eqnarray} +// A(q)y(n) = B(q)u(n-k) + e(t) +// end{eqnarray} +// </latex> +// IV4 model is SISO type model. It is unaffected by color of the noise. Four steps used in IV4 model design. First step is the generation of the ARX model. +// Second step uses the ARX model to generate the instrument variable matrix.Next steps uses the residual to generate a higher order model coefficient. +// In final step uses the AR model coefficient to filter the input and output data and feed it to the IV model. +// sys ,an idpoly type class, have different fields that contains estimated coefficients, sampling time, time unit and other estimated data in Report object. +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.2];b = [0 0.2 0.3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// ioData = iddata(y,u,0.1) +// sys = iv4(ioData,[2,2,1]) +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.2];b = [0 0.2 0.3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// ioData = [y,u] +// sys = iv4(ioData,[2,2,1]) +// +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + [lhs, rhs] = argn(0) - plantData = varargin(1) orderData = varargin(2) na = orderData(1);nb = orderData(2) - // arranging na ,nb,nk + // arranging na ,nb,nk if size(orderData,"*") == 2 then nk = 1 elseif size(orderData,'*') == 3 then @@ -12,7 +55,7 @@ function varargout = iv4(varargin) end nb1 = nb + nk - 1 n = max(na, nb1) - // arranging the plant data + // arranging the plant data if typeof(plantData) == 'constant' then Ts = 1;unitData = 'second' elseif typeof(plantData) == 'iddata' then @@ -20,9 +63,9 @@ function varargout = iv4(varargin) plantData = [plantData.OutputData plantData.InputData] end noOfSample = size(plantData,'r') - // finding the iv model + // finding the iv model ivTest = iv(plantData,[na nb nk]); - // residual + // residual [aTemp,bTemp,cTemp] = pe(plantData,ivTest); Lhat = ar(aTemp,na+nb); x = sim(plantData(:,2),ivTest); @@ -30,17 +73,17 @@ function varargout = iv4(varargin) Yf = filter(Lhat.a,Lhat.b,[plantData(:,1);zeros(n,1)]); phif = zeros(noOfSample,na+nb) psif = zeros(noOfSample,na+nb) - // arranging samples of y matrix + // arranging samples of y matrix for ii = 1:na phif(ii+1:ii+noOfSample,ii) = -yData psif(ii+1:ii+noOfSample,ii) = -x end - // arranging samples of u matrix + // arranging samples of u matrix for ii = 1:nb phif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData psif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData end - // passing it through the filters + // passing it through the filters for ii = 1:na+nb phif(:,ii) = filter(Lhat.a,Lhat.b,phif(:,ii)); psif(:,ii) = filter(Lhat.a,Lhat.b,psif(:,ii)); @@ -55,12 +98,12 @@ function varargout = iv4(varargin) vcov = sigma2 * pinv((phif)' * phif) t = idpoly([1; theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,Ts) - // estimating the other parameters - [temp1,temp2,temp3] = predict(z,t) - [temp11,temp22,temp33] = pe(z,t) + // estimating the other parameters + [temp1,temp2,temp3] = predict(plantData,t) + [temp11,temp22,temp33] = pe(plantData,t) estData = calModelPara(temp1,temp11,na+nb) - //pause + // pause t.Report.Fit.MSE = estData.MSE t.Report.Fit.FPE = estData.FPE t.Report.Fit.FitPer = estData.FitPer @@ -69,7 +112,7 @@ function varargout = iv4(varargin) t.Report.Fit.nAIC = estData.nAIC t.Report.Fit.BIC = estData.BIC t.TimeUnit = unitData - //sys = t + // sys = t varargout(1) = t - //varargout(1) = idpoly([1; -theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,Ts) + // varargout(1) = idpoly([1; -theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,Ts) endfunction diff --git a/misdata.sci b/misdata.sci index c1ced50..5bd6340 100644 --- a/misdata.sci +++ b/misdata.sci @@ -1,11 +1,37 @@ function varargout = misdata(varargin) + +// Recover Missing Data by Interpolation +// +// Calling Sequence +// data = misdata(plantData) +// Parameters +// plantData : iddata type object with missing data as Nan +// data : iddata type object with interpolated data +// +// Description +// misdata function recovers the experimental missing plant time series data by linear interpolation. +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.2];b = [0 0.2 0.3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// u(100:105) = %nan;u(300:100:1000) = %nan +// y(420:422) = %nan;y(555:100:1024) = %nan +// plantData = iddata(y,u,0.1) +// data = misdata(plantData) +// +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + + [lhs,rhs] = argn(0) -//------------------------------------------------------------------------------ -// checking the number of inputs +// ------------------------------------------------------------------------------ +// checking the number of inputs if rhs <> 1 then error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"misdata")) end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ ioData = varargin(1) if typeof(ioData) <> "iddata" then error(msprintf(gettext("%s: Plant input data must be ""iddata"" type. "),"misdata")) @@ -16,14 +42,14 @@ function varargout = misdata(varargin) endfunction function varargout = linearINTRP(matData,Ts) - // looking for overall nan values + // looking for overall nan values nanData = isnan(matData);nanIndex = find(nanData == %T) if ~size(nanIndex,'*') then varargout(1) = matData else tempMat = [] matSize = size(matData,'r') - // looking for nan in each column + // looking for nan in each column for ii = 1:size(matData,'c') nanData = isnan(matData(:,ii));nanIndex = find(nanData == %T); if ~size(nanData,'*') then diff --git a/nInputData.sci b/nInputData.sci index 59a3b64..8d5a481 100644 --- a/nInputData.sci +++ b/nInputData.sci @@ -1,15 +1,50 @@ function varargout = nInputSeries(varargin) - [lhs rhs] = argn(0) +// Count the number of input output data series +// +// Calling Sequence +// ioData = nInputSeries(plantData) +// +// Parameters +// plantData : iddata type module +// ioData : structure type variable having input and output object +// +// Description +// plantData must be iddata type. The function returns the number of input and output data series in struct form. It has two fields name as input and output. +// +// Examples +// plantData = iddata(rand(100,3),rand(100,2),0.1) +// ioData = nInputSeries(plantData) +// +// Examples +// plantData = iddata(rand(100,2),[],0.1) +// ioData = nInputSeries(plantData) +// +// Examples +// plantData = iddata([],rand(100,3),0.1) +// ioData = nInputSeries(plantData) +// +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar +// + [lhs rhs] = argn(0) if rhs <> 1 then - error(msprintf(gettext("%s: Wrong number of input arguments."),"iddataplot")) + error(msprintf(gettext("%s: Wrong number of input arguments."),"nInputSeries")) end iddataData = varargin(1) if typeof(iddataData) <> 'iddata' then error(msprintf(gettext("%s:Wrong type for input argument %d: ""iddata"" expected."),"nInputSeries",1)) end + t = struct('input',0,'output',0) + if ~size(iddataData.InputData,'*') then - varargout(1) = 1 + t.input = 0 else - varargout(1) = size(iddataData.InputData,'c') + t.input = size(iddataData.InputData,'c') end + if ~size(iddataData.OutputData,'*') then + t.output = 0 + else + t.output = size(iddataData.InputData,'c') + end + varargout(1) = t endfunction @@ -1,19 +1,49 @@ -// Copyright (C) 2015 - IIT Bombay - FOSSEE -// -// This file must be used under the terms of the CeCILL. -// This source file is licensed as described in the file COPYING, which -// you should have received as part of this distribution. The terms -// are also available at -// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt -// Authors: Harpreet, Ashutosh -// Organization: FOSSEE, IIT Bombay function sys = oe(varargin) -// Estimates Discrete time BJ 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 +// Parameters Estimation of OE(Output-Error) model using Input Output time-domain data +// +// Calling Sequence +// sys = oe(ioData,[nb nf nk]) +// +// Parameters +// ioData : iddata or [outputData inputData] ,matrix of nx2 dimensions, type plant data +// nb : non-negative integer number specified as order of the polynomial A(z^-1) +// nf : non-negative integer number specified as order of the polynomial B(z^-1)+1 +// nk : non-negative integer number specified as input output delay, Default value is 1 +// sys : idpoly type polynomial have estimated coefficients of B(z^-1),f(z^-1) polynomials +// +// Description +// Fit OE model on given input output data +// The mathematical equation of the OE model +// <latex> +// begin{eqnarray} +// y(n) = \frac {B(q)}{D(q)}u(n) + e(t) +// end{eqnarray} +// </latex> +// It is SISO type model. It minimizes the sum of the squares of nonlinear functions using Levenberg-Marquardt algorithm. +// sys ,an idpoly type class, have different fields that contains estimated coefficients, sampling time, time unit and other estimated data in Report object. +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 2];b = [0 2 3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// plantData = iddata(y,u,0.1) +// sys = oe(plantData,[2,2,1]) +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 2];b = [0 2 3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// plantData = [y,u] +// sys = oe(plantData,[2,2,1]) +// +// Authors +// Ashutosh Kumar Bhargava, Harpreet,Inderpreet + + [lhs , rhs] = argn(); if ( rhs < 2 ) then @@ -38,7 +68,7 @@ function sys = oe(varargin) errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "oe_2"); error(errmsg); end -// +// n = varargin(2) if (size(n,"*")<2| size(n,"*")>3) then errmsg = msprintf(gettext("%s: The order and delay matrix [nb nf nk] should be of size [2 4]"), "oe_2"); @@ -49,9 +79,9 @@ function sys = oe(varargin) errmsg = msprintf(gettext("%s: values of order and delay matrix [nb nf nk] should be nonnegative integer number "), "oe_2"); error(errmsg); end -// +// nb= n(1); nf = n(2); -// +// if (size(n,"*") == 2) then nk = 1 else @@ -80,7 +110,7 @@ function sys = oe(varargin) [temp11,temp22,temp33] = pe(z,t) estData = calModelPara(temp1,temp11,n(1)+n(2)) - //pause + // pause t.Report.Fit.MSE = estData.MSE t.Report.Fit.FPE = estData.FPE t.Report.Fit.FitPer = estData.FitPer @@ -90,8 +120,8 @@ function sys = oe(varargin) t.Report.Fit.BIC = estData.BIC t.TimeUnit = unit sys = t - //sys = t - //sys.TimeUnit = unit + // sys = t + // sys.TimeUnit = unit endfunction function yhat = _objoefun(UDATA,YDATA,x,nf,nb,nk) @@ -1,8 +1,38 @@ function varargout = pe(varargin) + +// K-steps ahead prediction error +// +// Calling Sequence +// pe(plantData,sys) +// pe(plantData,sys,k) +// [yData,tData,fData] = pe(plantData,sys) +// [yData,tData,fData] = pe(plantData,sys,k) +// Parameters +// plantData : iddata type or nx2 matrix +// sys : idpoly type polynomial +// k : non-neagtive integer prediction step,default value is 1 +// yData : k step ahead prediction error +// tData : time series data +// fData : initial state +// Description +// pe function estimate the differencr between k step ahead predicted and the actual output response. Initial conditions of the dynamic model is zero. +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// y = sim(u,sys)+rand(1024,1) +// plantData = iddata(y,u,0.1) +// pe(plantData,sys) +// figure();clf(); +// k = 5 +// pe(plantData,sys,k) +// Authors +// Ashutosh Kumar Bhargava + [lhs,rhs] = argn(0) -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ data = varargin(1) model = varargin(2) if rhs == 3 then @@ -10,7 +40,7 @@ function varargout = pe(varargin) elseif rhs == 2 then kStep = 1 end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ // k step analysis if typeof(kStep) <> 'constant' || isnan(kStep) then error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"pe")) @@ -23,15 +53,15 @@ function varargout = pe(varargin) if size(kStep,'*') <> 1 || (ceil(kStep)-kStep) then error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"pe")) end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ // checking the plant model if typeof(model) ~= 'idpoly' then error(msprintf(gettext("%s:Plant model must be ""idpoly"" type.\n"),"pe")) end modelSampleTime = model.Ts modelTimeUnit = model.TimeUnit -//------------------------------------------------------------------------------ -//checking the data type +// ------------------------------------------------------------------------------ +// checking the data type if typeof(data) <> 'iddata' && typeof(data) <> 'constant' then error(msprintf(gettext("%s:Sample data must be ""iddata"" type or ""n x 2"" matrix type.\n"),"pe")) end @@ -43,7 +73,7 @@ function varargout = pe(varargin) plantSampleTime = data.Ts plantTimeUnit = data.TimeUnit data = [data.OutputData data.InputData] - //disp('iddata') + // disp('iddata') elseif typeof(data) == 'constant' then if size(data,'c') ~= 2 then error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"pe")) @@ -51,7 +81,7 @@ function varargout = pe(varargin) plantSampleTime = model.Ts plantTimeUnit = model.TimeUnit end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ // comparing the sampling time if modelSampleTime-plantSampleTime <> 0 then error(msprintf(gettext("%s:The sample time of the model and plant data must be equal.\n"),"pe")) @@ -61,7 +91,7 @@ function varargout = pe(varargin) else error(msprintf(gettext("%s:Time unit of the model and plant data must be equal.\n"),"pe")) end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ sampleLength = size(data,'r') // one step ahead prediction [y1 ,x0] = predict(data,model) @@ -89,7 +119,7 @@ function varargout = pe(varargin) timeData = (modelSampleTime:modelSampleTime:(sampleLength)*modelSampleTime)' pseudoData = size(eCapData,'r') eCapData = eCapData(abs(pseudoData-sampleLength)+1:$) - //pause + // pause if lhs == 1 then clf() plot(timeData,eCapData) @@ -98,7 +128,7 @@ function varargout = pe(varargin) tempTimeUnit = 'Time('+modelTimeUnit+')' xtitle('Predicted Response',tempTimeUnit,'y') xgrid - //pause + // pause varargout(1) = 0 elseif lhs == 2 then diff --git a/predict.sci b/predict.sci index fcd4b23..2eaba3c 100644 --- a/predict.sci +++ b/predict.sci @@ -1,17 +1,46 @@ -//[y,x0] = predict(data,idpoly,k) -//References -//Digital Control(11.1.2) by Kanna M.Moudgalya -//System Identification Theory for User Second Edition (3.2) by Lennart Ljung -// Code Aurthor - Ashutosh Kumar Bhargava function varargout = predict(varargin) + +// K-steps ahead output predictor +// +// Calling Sequence +// predict(plantData,sys) +// predict(plantData,sys,k) +// [yData,tData,fData] = predict(plantData,sys) +// [yData,tData,fData] = predict(plantData,sys,k) +// +// Parameters +// plantData : iddata type or nx2 matrix +// sys : idpoly type polynomial +// k : non-neagtive integer prediction step +// yData : k step ahead predicted output response,default value is 1 +// tData : time series data +// fData : initial state +// +// Description +// predict function consider the inital conditions as zero and predict the k step ahead output response of the sys ,idpoly type, model. +// +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// y = sim(u,sys)+rand(1024,1) +// plantData = iddata(y,u,0.1) +// predict(plantData,sys) +// figure();clf(); +// k = 5 +// predict(plantData,sys,k) +// +// Authors +// Ashutosh Kumar Bhargava + [lhs,rhs] = argn(0) -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ // checking the number of inputs if rhs < 2 || rhs > 3 then error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"predict")) end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ data = varargin(1) model = varargin(2) if rhs == 3 then @@ -20,7 +49,7 @@ function varargout = predict(varargin) kStep = 1 end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ // k step analysis if typeof(kStep) <> 'constant' || isnan(kStep) then @@ -34,15 +63,15 @@ function varargout = predict(varargin) if size(kStep,'*') <> 1 || (ceil(kStep)-kStep) then error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"predict")) end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ // checking the plant model if typeof(model) ~= 'idpoly' then error(msprintf(gettext("%s:Plant model must be ""idpoly"" type.\n"),"predict")) end modelSampleTime = model.Ts modelTimeUnit = model.TimeUnit -//------------------------------------------------------------------------------ -//checking the data type +// ------------------------------------------------------------------------------ +// checking the data type if typeof(data) <> 'iddata' && typeof(data) <> 'constant' then error(msprintf(gettext("%s:Sample data must be ""iddata"" type or ""n x 2"" matrix type.\n"),"predict")) end @@ -54,7 +83,7 @@ function varargout = predict(varargin) plantSampleTime = data.Ts plantTimeUnit = data.TimeUnit data = [data.OutputData data.InputData] - //disp('iddata') + // disp('iddata') elseif typeof(data) == 'constant' then if size(data,'c') ~= 2 then error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"predict")) @@ -62,7 +91,7 @@ function varargout = predict(varargin) plantSampleTime = model.Ts plantTimeUnit = model.TimeUnit end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ // comparing the sampling time if modelSampleTime-plantSampleTime <> 0 then error(msprintf(gettext("%s:The sample time of the model and plant data must be equal.\n"),"predict")) @@ -72,14 +101,14 @@ function varargout = predict(varargin) else error(msprintf(gettext("%s:Time unit of the model and plant data must be equal.\n"),"predict")) end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ // ckecking the k step size. if it greater than number of sample size then the // k step will become 1 if kStep >= size(data,'r') then kStep = 1 end -//------------------------------------------------------------------------------ - //storing the plant data +// ------------------------------------------------------------------------------ + // storing the plant data // B(z) C(z) // y(n) = ---------- u(n) + ---------- e(n) // A(z)*F(z) A(z)*D(z) @@ -114,17 +143,17 @@ function varargout = predict(varargin) Wkq1.num = Wkq1.num/tempWkq1(1) Wkq1.den = Wkq1.den/tempWkq1(1) end - //pause -//------------------------------------------------------------------------------ + // pause +// ------------------------------------------------------------------------------ // storing the plant data uCoeff = coeff(WkqGq.num*Wkq1.den) yCoeff = coeff(WkqGq.den*Wkq1.num) yCapCoeff = coeff(WkqGq.den*Wkq1.den) - //pause + // pause lengthuCoeff = length(uCoeff) lengthyCoeff = length(yCoeff) lengthyCapCoeff = length(yCapCoeff) -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ // keeping initial conditions equal to zero uData = zeros(lengthuCoeff,1) yData = zeros(lengthyCoeff,1) @@ -132,7 +161,7 @@ function varargout = predict(varargin) uData = [uData;data(:,2)] yData = [yData;data(:,1)] sampleData = size(data,'r') - //pause + // pause // reversing the coefficients if ~size(uCoeff,'*') then uCoeff = 0 @@ -149,9 +178,9 @@ function varargout = predict(varargin) else yCapCoeff = -yCapCoeff(lengthyCapCoeff:-1:2) end - //pause + // pause for ii = 1:sampleData+1 - //pause + // pause if ~size(uData(ii:ii+lengthuCoeff-1),'*') then tempu = 0 else @@ -1,6 +1,42 @@ function varargout = rarx(varargin) + +// Parameters Estimation of ARX model by recursive method +// +// Calling Sequence +// sys = rarx(ioData,[na nb nk],lambda) +// Parameters +// ioData : iddata or [outputData inputData] ,matrix of nx2 dimensions, type plant data +// na : non-negative integer number specified as order of the polynomial A(z^-1) +// nb : non-negative integer number specified as order of the polynomial B(z^-1)+1 +// nk : non-negative integer number specified as input output delay, Default value is 1 +// lambda : Forgetting factor,Default value is 0.95 +// sys : idpoly type polynomial have estimated coefficients of A(z^-1) and B(z^-1) polynomials +// +// Description +// Fit RARX model on given input output data +// RARX model is SISO type model. It uses recursive weighted least-squares algorithm to estimate the coefficient of ARX model +// sys is a struct type variable output contains data about theta and yhat. +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.2];b = [0 0.2 0.3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// ioData = iddata(y,u,0.1) +// sys = rarx(ioData,[2,2,1]) +// +// Examples +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.2];b = [0 0.2 0.3]; +// model = idpoly(a,b,'Ts',0.1) +// y = sim(u,model) + rand(length(u),1) +// ioData = [y,u] +// sys = rarx(ioData,[2,2,1]) +// +// Authors +// Ashutosh Kumar Bhargava, Bhushan Manjarekar + [lhs,rhs] = argn(0) -// plantData = varargin(1) orderData = varargin(2) na = orderData(1);nb = orderData(2) @@ -44,7 +80,7 @@ function varargout = rarx(varargin) for ii = 1:nb tempData(ii+nk:ii+N+nk-1,ii+na) = plantData(:,2) end - //tempData = [zeros(1,na+nb);tempData] + // tempData = [zeros(1,na+nb);tempData] tempData = tempData(1:N+1,:) for ii = 1:N diff --git a/resid.sci b/resid.sci new file mode 100644 index 0000000..1184c8d --- /dev/null +++ b/resid.sci @@ -0,0 +1,86 @@ +function resid(varargin) + +// plot residual analysis +// +// Calling Sequence +// resid(plantData,sys) +// Parameters +// plantData : iddata type or nx2 matrix +// sys : idpoly type polynomial +// Description +// resid function compute the residuals of the model with its plant data by taking auto-correlation of the one step ahead prediction error and cross-correlations +// in between with one step ahead prediction error with plant input. +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// y = sim(u,sys)+rand(1024,1) +// plantData = iddata(y,u,0.1) +// resid(plantData,sys) +// Authors +// Ashutosh Kumar Bhargava + + [lhs rhs] = argn(0) + // Checking the number of inputs + if rhs <> 2 then + error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"resid")) + end + // Checking the type of first argument + if typeof(varargin(1)) <> 'constant' && typeof(varargin(1)) <> 'iddata' then + error(msprintf(gettext("%s:Incompatible model data.\n"),"resid")) + end + // checking the type of second argument + if typeof(varargin(2)) <> 'idpoly' then + error(msprintf(gettext("%s:Wrong type for second argument,""idpoly"" type expected.\n"),"resid")) + end + // checking the input dimensions + yData = [] + uData = [] + tempData = varargin(1) + if typeof(varargin(1)) == 'iddata' then + yData = tempData.OutputData + uData = tempData.InputData + if ~size(yData,'*') || ~size(uData,'*') then + error(msprintf(gettext("%s:Model input and output numbers must be equal to the plant data.\n"),"resid")) + end + elseif typeof(varargin(1)) == 'constant'then + if isrow(tempData) || iscolumn(tempData) || length(size(tempData)) <> 2 || size(tempData,'c') <> 2 then + error(msprintf(gettext("%s:Model input and output numbers must be equal to the plant data.\n"),"resid")) + end + yData = tempData(:,1) + uData = tempData(:,2) + tempData = iddata(yData,uData,varargin(2).Ts) + end + + // [a,b,c] = compare(tempData,varargin(2)) + [a,x0] = pe(tempData,varargin(2)) +// ------------------------------------------------------------------------------ + errorData = a + eAutoCorr = xcorr(errorData,25,'coeff')// 'biased') + uAutoCorr = xcorr(uData,25,'coeff')// 'biased') + eAutoCorrMax = max(eAutoCorr) + eAutoCorr = eAutoCorr/eAutoCorrMax +// ------------------------------------------------------------------------------ + + // confidence interval + coi1 = 2.58/sqrt(length(yData)) + rect = [0;coi1;25;2*coi1]; + // disp(rect) + subplot(2,1,1); + xrects(rect,7) + stem([0:25]',eAutoCorr(26:length(eAutoCorr)),5) + xrect(rect) + coi2 = 2.58*sqrt(eAutoCorrMax+2*(eAutoCorr(27:50)'*uAutoCorr(27:50)*eAutoCorrMax))/sqrt(length(yData)) + rect = [-25;coi2;50;2*coi2]; + euCrossCorr = xcorr(errorData,uData,25,'biased') + // pause + xtitle('Autocorrelation of Residuals','lag') + subplot(2,1,2); + xrects(rect,7) + // disp(rect) + stem([-25:25]',euCrossCorr,5) + xrect(rect) + xtitle('Crosscorrelation of Residuals and Input','lag') + // figure + // plot2d3(euCrossCorr) +endfunction diff --git a/scilab_error.sci b/scilab_error.sci new file mode 100644 index 0000000..9eadde0 --- /dev/null +++ b/scilab_error.sci @@ -0,0 +1,20 @@ +// ==================================================================== +// Template toolbox_skeleton +// This file is released under the 3-clause BSD license. See COPYING-BSD. +// ==================================================================== +// +// +function scilab_error(varargin) + + argSize = size(varargin); + + //in toolboxes, use "_d" or "dgettext" to your localized messages + if argSize <> 1 then + error(999, msprintf(_d("toolbox_skeleton", "%s: I''m waiting for only one argument.\n"), "scilab_error")); + end + + if argSize == 1 then + error(999, msprintf(dgettext("toolbox_skeleton", "%s: Yeah! %d is a good number of arguments but I prefer fail, sorry.\n"), "scilab_error", 1)); + end +endfunction +// ==================================================================== diff --git a/scilab_sum.sci b/scilab_sum.sci new file mode 100644 index 0000000..b9c5fa0 --- /dev/null +++ b/scilab_sum.sci @@ -0,0 +1,10 @@ +// ==================================================================== +// Template toolbox_skeleton +// This file is released under the 3-clause BSD license. See COPYING-BSD. +// ==================================================================== +// +// +function s = scilab_sum(valA,valB) + s = valA + valB; +endfunction +// ==================================================================== @@ -1,12 +1,32 @@ function varargout = sim(varargin) + +// Simulate response of idpoly type polynomials +// +// Calling Sequence +// yData = sim(uData,sys) +// Parameters +// sys : idpoly type polynomial +// uData : plant input data in iddata class or nx1 matrix +// yData : simulated plant output data nx1 matrix +// Description +// sim function returnes the simulated response of idpoly type dynamic function for the given input data +// +// Examples +// uData = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// yData = sim(uData,sys) +// Authors +// Ashutosh Kumar Bhargava + [lhs,rhs] = argn(0) - //checking the number of inputs + // checking the number of inputs if rhs < 2 || rhs > 3 then error(msprintf(gettext("%s: Unexpected number of input arguments "),"sim")) end modelData = varargin(2) inputData = varargin(1) - //checking the first input + // checking the first input if typeof(inputData) <> "constant" && typeof(inputData) <> "iddata" then error(msprintf(gettext("%s: Plant input data must be ""iddata"" type or ""double vector"". "),"sim")) end @@ -16,11 +36,11 @@ function varargout = sim(varargin) if ~iscolumn(inputData) then error(msprintf(gettext("%s: Plant input data must be ""double column vector"". "),"sim")) end - //checking the plant model type + // checking the plant model type if typeof(modelData) <> "idpoly" then error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"sim")) end - //adding noise + // adding noise if rhs == 3 then noiseFlag = varargin(3) if typeof(noiseFlag) <> 'boolean' then @@ -29,7 +49,7 @@ function varargout = sim(varargin) else noiseFlag = %F end - //adding noise of zero mean and 1 as standard deviation + // adding noise of zero mean and 1 as standard deviation if noiseFlag then numberOfSamples = size(inputData,'r') R = = grand(numberOfSamples,1,"nor",0,1) @@ -1,4 +1,26 @@ function varargout = spa(varargin) +// Use spectral analysis to estimate frequency response +// +// Calling Sequence +// frdData = spa(plantData,winSize,Freq) +// Parameters +// plantData : iddata type +// winSize : non-neative integer number +// Freq : frequency points to evaluate the response +// frdData : frd type object +// Description +// spa function does the estimation of the frequency response of iddata of a plant using the spectral analysis. Hanning window is used in spectral +// analysis. Default size of the window is minimum of the number of sample divided by 10 or 30. The default frquency point is %pi x (1:128)/(sampling time x 128). +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// u = idinput(1024,'PRBS',[0 1/20],[-1 1]) +// y = sim(u,sys)+rand(1024,1) +// plantData = iddata(y,u,0.1) +// frdData = spa(plantData) +// Authors +// Bhushan Manjarekar, Ashutosh Kumar Bhargava + [lhs , rhs] = argn(); if ( rhs < 1 || rhs > 3 ) then errmsg = msprintf(gettext("%s: Wrong number of input arguments" ),"spa"); @@ -10,7 +32,7 @@ function varargout = spa(varargin) plantData = varargin(1) windowSize = %F inputFreq = %F -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ // arranging the plant data inputData = plantData.InputData; if ~size(inputData,"*") then @@ -20,7 +42,7 @@ function varargout = spa(varargin) if ~size(outputData,"*") then error(msprintf(gettext("%s:Output data must be non-zero vector. \n"),"spa")) end -//------------------------------------------------------------------------------ +// ------------------------------------------------------------------------------ N = size(inputData,'r') nout = size(outputData,'c') nin = size(inputData,'c') @@ -43,10 +65,10 @@ function varargout = spa(varargin) phiY = spaCalculation(inputFreq,Ryy,M) temp = phiY for jj = 1:nin - phiYU = spaCalculation(inputFreq,Ryu,M)//sapply(freq, cov2spec, Ryu[i, j, ], M) + phiYU = spaCalculation(inputFreq,Ryu,M)// sapply(freq, cov2spec, Ryu[i, j, ], M) phiU = spaCalculation(inputFreq,Ruu,M) G = [G phiYU./phiU] - //pause + // pause temp = temp - phiYU .* conj(phiYU)./phiU end spect = [spect temp] diff --git a/stepest.sci b/stepest.sci index 6310971..3c50766 100644 --- a/stepest.sci +++ b/stepest.sci @@ -1,6 +1,31 @@ function varargout = stepest(varargin) + +// Estimate step response and plot of idpoly type model +// +// Calling Sequence +// stepest(sys) +// stepData = stepest(sys,flag) +// Parameters +// sys : idpoly type polynomial +// flag : boolean type variable,default value is false(%F) +// stepData : stores step response if the flag value is true +// Description +// stepest function estimate and plot the impulse response of idpoly type function. +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// stepest(sys); +// +// Examples +// a = [1 0.2];b = [0 0.2 0.3]; +// sys = idpoly(a,b,'Ts',0.1) +// flag = %T +// stepData = stepest(sys,flag) +// Authors +// Ashutosh Kumar Bhargava + [lhs,rhs] = argn(0) - //checking the number of inputs + // checking the number of inputs if rhs > 2 then error(msprintf(gettext("%s: Unexpected number of input arguments "),"stepest")) end @@ -8,7 +33,7 @@ function varargout = stepest(varargin) if typeof(modelData) <> "idpoly" then error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"stepest")) end - //adding noise + // adding noise if rhs == 2 then noiseFlag = varargin(2) if typeof(noiseFlag) <> 'boolean' then |