added scilabs files

1 files changed, 192 insertions, 0 deletions

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