From e5e316e1958e27696d7670e2492992d34ff38b68 Mon Sep 17 00:00:00 2001 From: ttt Date: Mon, 9 Jul 2018 16:54:44 +0530 Subject: added scilabs files --- predict.sci | 192 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 192 insertions(+) create mode 100644 predict.sci (limited to 'predict.sci') diff --git a/predict.sci b/predict.sci new file mode 100644 index 0000000..fcd4b23 --- /dev/null +++ b/predict.sci @@ -0,0 +1,192 @@ +//[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) + [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 + kStep = varargin(3) + 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"),"predict")) + 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"),"predict")) + 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"),"predict")) + 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"),"predict")) + 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"),"predict")) + end + 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")) + 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"),"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 + // 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) + + 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 + end + 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) + end + //pause +//------------------------------------------------------------------------------ + // storing the plant data + uCoeff = coeff(WkqGq.num*Wkq1.den) + yCoeff = coeff(WkqGq.den*Wkq1.num) + yCapCoeff = coeff(WkqGq.den*Wkq1.den) + //pause + lengthuCoeff = length(uCoeff) + lengthyCoeff = length(yCoeff) + lengthyCapCoeff = length(yCapCoeff) +//------------------------------------------------------------------------------ + // keeping initial conditions equal to zero + uData = zeros(lengthuCoeff,1) + yData = zeros(lengthyCoeff,1) + yCapData = zeros(lengthyCapCoeff-1,1) + uData = [uData;data(:,2)] + yData = [yData;data(:,1)] + sampleData = size(data,'r') + //pause + // reversing the coefficients + if ~size(uCoeff,'*') then + uCoeff = 0 + else + uCoeff = uCoeff(lengthuCoeff:-1:1) + end + if ~size(yCoeff) then + yCoeff = 0 + else + yCoeff = yCoeff(lengthyCoeff:-1:1) + end + if ~size(yCapCoeff,'*') then + yCapCoeff = 0 + else + yCapCoeff = -yCapCoeff(lengthyCapCoeff:-1:2) + end + //pause + for ii = 1:sampleData+1 + //pause + if ~size(uData(ii:ii+lengthuCoeff-1),'*') then + tempu = 0 + else + tempu = uCoeff*uData(ii:ii+lengthuCoeff-1); + end + if ~size(yData(ii:ii+lengthyCoeff-1),'*') + tempy = 0 + else + tempy = yCoeff*yData(ii:ii+lengthyCoeff-1); + end + if ~size(yCapData(ii:ii+lengthyCapCoeff-2),'*') then + tempyCap = 0 + else + tempyCap = yCapCoeff*yCapData(ii:ii+lengthyCapCoeff-2); + end + yCapData = [yCapData;tempu+tempy+tempyCap]; + end + // pause + extraSample = abs(size(yCapData,'r')-sampleData) + yCapData = yCapData(extraSample+1:$) + timeData = ((modelSampleTime:modelSampleTime:(size(yCapData,'r')*modelSampleTime))'); + if lhs == 1 then + clf() + plot(timeData,yCapData) + axisData = gca() + tempTimeUnit = 'Time('+modelTimeUnit+')' + xtitle('Predicted Response',tempTimeUnit,'y') + xgrid + varargout(1) = 0 + elseif lhs == 2 then + varargout(1) = yCapData + varargout(2) = 0 + elseif lhs == 3 then + varargout(1) = yCapData + varargout(2) = timeData + varargout(3) = 0 + end +endfunction -- cgit