diff options
Diffstat (limited to 'compare.sci')
| -rw-r--r-- | compare.sci | 250 |
1 files changed, 124 insertions, 126 deletions
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 + |