diff options
Diffstat (limited to 'iv4.sci')
| -rw-r--r-- | iv4.sci | 75 |
1 files changed, 75 insertions, 0 deletions
@@ -0,0 +1,75 @@ +function varargout = iv4(varargin) + [lhs, rhs] = argn(0) + + plantData = varargin(1) + orderData = varargin(2) + na = orderData(1);nb = orderData(2) + // arranging na ,nb,nk + if size(orderData,"*") == 2 then + nk = 1 + elseif size(orderData,'*') == 3 then + nk = orderData(3) + end + nb1 = nb + nk - 1 + n = max(na, nb1) + // arranging the plant data + if typeof(plantData) == 'constant' then + Ts = 1;unitData = 'second' + elseif typeof(plantData) == 'iddata' then + Ts = plantData.Ts;unitData = plantData.TimeUnit + plantData = [plantData.OutputData plantData.InputData] + end + noOfSample = size(plantData,'r') + // finding the iv model + ivTest = iv(plantData,[na nb nk]); + // residual + [aTemp,bTemp,cTemp] = pe(plantData,ivTest); + Lhat = ar(aTemp,na+nb); + x = sim(plantData(:,2),ivTest); + yData = plantData(:,1);uData = plantData(:,2) + 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 + for ii = 1:na + phif(ii+1:ii+noOfSample,ii) = -yData + psif(ii+1:ii+noOfSample,ii) = -x + end + // 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 + for ii = 1:na+nb + phif(:,ii) = filter(Lhat.a,Lhat.b,phif(:,ii)); + psif(:,ii) = filter(Lhat.a,Lhat.b,psif(:,ii)); + end + lhs = psif'*phif + lhsinv = pinv(lhs) + theta = lhsinv * (psif)' * Yf + ypred = (phif * theta) + ypred = ypred(1:size(yData,'r')) + e = yData - ypred + sigma2 = norm(e)^2/(size(yData,'r') - na - nb) + 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) + + estData = calModelPara(temp1,temp11,na+nb) + //pause + t.Report.Fit.MSE = estData.MSE + t.Report.Fit.FPE = estData.FPE + t.Report.Fit.FitPer = estData.FitPer + t.Report.Fit.AIC = estData.AIC + t.Report.Fit.AICc = estData.AICc + t.Report.Fit.nAIC = estData.nAIC + t.Report.Fit.BIC = estData.BIC + t.TimeUnit = unitData + //sys = t + varargout(1) = t + //varargout(1) = idpoly([1; -theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,Ts) +endfunction |