diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/data.R | 117 | ||||
| -rw-r--r-- | R/estUtil.R | 232 | ||||
| -rw-r--r-- | R/estpoly.R | 636 | ||||
| -rw-r--r-- | R/idframe.R | 280 | ||||
| -rw-r--r-- | R/idinput.R | 156 | ||||
| -rw-r--r-- | R/imports.R | 5 | ||||
| -rw-r--r-- | R/ioNamesData.R | 129 | ||||
| -rw-r--r-- | R/iv.R | 157 | ||||
| -rw-r--r-- | R/nonparam.R | 271 | ||||
| -rw-r--r-- | R/poly.R | 196 | ||||
| -rw-r--r-- | R/predict.R | 120 | ||||
| -rw-r--r-- | R/preprocess.R | 192 | ||||
| -rw-r--r-- | R/rarx.R | 79 | ||||
| -rw-r--r-- | R/readData.R | 73 | ||||
| -rw-r--r-- | R/sim.R | 87 | ||||
| -rw-r--r-- | R/sysid.R | 0 | ||||
| -rw-r--r-- | R/util.R | 53 |
17 files changed, 2783 insertions, 0 deletions
diff --git a/R/data.R b/R/data.R new file mode 100644 index 0000000..89b65ea --- /dev/null +++ b/R/data.R @@ -0,0 +1,117 @@ +#' Data simulated from an ARMAX model +#' +#' This dataset contains 2555 samples simulated from the following ARMAX model: +#' \deqn{ +#' y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + +#' \frac{1-0.3q^{-1}}{1 - 0.5q^{-1}} e[k] +#' } +#' +#' The model is simulated with a 2555 samples long full-band PRBS input. +#' The noise variance is set to 0.1 +#' +#' @format an \code{idframe} object with 2555 samples, one input and one +#' output +#' +"armaxsim" + +#' Data simulated from an ARX model +#' +#' This dataset contains 2555 samples simulated from the following ARX model: +#' \deqn{ +#' y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + +#' \frac{1}{1 - 0.5q^{-1}} e[k] +#' } +#' +#' The model is simulated with a 2555 samples long full-band PRBS input. +#' The noise variance is set to 0.1 +#' +#' @format an \code{idframe} object with 2555 samples, one input and one +#' output +#' +"arxsim" + +#' Data simulated from an BJ model +#' +#' This dataset contains 2046 samples simulated from the following BJ model: +#' \deqn{ +#' y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + +#' \frac{1+0.2q^{-1}}{1 - 0.3q^{-1}} e[k] +#' } +#' +#' The model is simulated with a 2046 samples long full-band PRBS input. +#' The noise variance is set to 0.1 +#' +#' @format an \code{idframe} object with 2046 samples, one input and one +#' output +#' +"bjsim" + +#' Continuous stirred tank reactor data (idframe) +#' +#' The Process is a model of a Continuous Stirring Tank Reactor, +#' where the reaction is exothermic and the concentration is +#' controlled by regulating the coolant flow. +#' \cr +#' +#' Inputs: q, Coolant Flow l/min +#' Outputs: +#' \describe{ +#' \item{Ca}{Concentration mol/l} +#' \item{T}{Temperature Kelvin}} +#' +#' @format an \code{idframe} object with 7500 samples, one input and two +#' outputs +"cstr" + +#' Continuous stirred tank reactor data (data.frame) +#' +#' The Process is a model of a Continuous Stirring Tank Reactor, +#' where the reaction is exothermic and the concentration is +#' controlled by regulating the coolant flow. +#' \cr +#' +#' Inputs: q, Coolant Flow l/min +#' Outputs: +#' \describe{ +#' \item{Ca}{Concentration mol/l} +#' \item{T}{Temperature Kelvin}} +#' +#' @format an \code{data.frame} object with 7500 rows and three columns: +#' q, Ca and T +#' +#' @source \url{ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/cstr.dat.gz} +"cstrData" + +#' Continuous stirred tank reactor data with missing values +#' +#' This dataset is derived from the \code{cstr} dataset with few samples +#' containing missing values, in one or all variables. It is used to +#' demonstrate the capabilities of the \code{misdata} routine. +#' +#' @format an \code{idframe} object with 7500 samples, one input and two +#' outputs +#' +#' @seealso \code{\link{cstr}}, \code{\link{misdata}} +"cstr_mis" + +#' Frequency response data +#' +#' This dataset contains frequency response data of an unknown SISO system. +#' +#' @format an \code{idfrd} object with response at 128 frequency points +"frd" + +#' Data simulated from an OE model +#' +#' This dataset contains 2555 samples simulated from the following OE model: +#' \deqn{ +#' y[k] = \frac{0.6q^{-2} - 0.2q^{-3}}{1 - 0.5q^{-1}} u[k] + e[k] +#' } +#' +#' The model is simulated with a 2555 samples long full-band PRBS input. +#' The noise variance is set to 0.1 +#' +#' @format an \code{idframe} object with 2555 samples, one input and one +#' output +#' +"oesim"
\ No newline at end of file diff --git a/R/estUtil.R b/R/estUtil.R new file mode 100644 index 0000000..084760c --- /dev/null +++ b/R/estUtil.R @@ -0,0 +1,232 @@ +# Implementation of the Levenberg Marquardt Algorithm +levbmqdt <- function(...,obj,theta0,N,opt){ + dots <- list(...) + + # Optimization Parameters + tol <- opt$tol; maxIter <- opt$maxIter + d <- opt$adv$LMinit; mu <- opt$adv$LMstep + + df <- N - dim(theta0)[1] + + # Initialize Algorithm + i <- 0 + l <- obj(theta=theta0,e=NULL,dots) + e <- l$e; grad <- l$grad + sumsq0 <- sum(e^2)/df + + # variable to count the number of times objective function is called + countObj <- 0 + sumSqDiff <- 9E-3*sumsq0 + + repeat{ + i=i+1 + + g <- 1/df*t(grad)%*%e + termPar <- norm(g,"2") + + repeat{ + # Update Parameters + H <- 1/df*t(grad)%*%grad + d*diag(dim(theta0)[1]) + Hinv <- solve(H); + + theta <- theta0 + Hinv%*%g + + # Evaulate sum square error + l <- obj(theta,e,dots) + sumsq <- sum(l$fn^2)/df + sumSqDiff <- sumsq0-sumsq + countObj <- countObj + 1 + + if(termPar < tol) break + # no major improvement + if(abs(sumSqDiff) < 0.01*sumsq0) break + + # If sum square error with the updated parameters is less than the + # previous one, the updated parameters become the current parameters + # and the damping coefficient is reduced by a factor of mu + if(sumSqDiff > 0){ + d <- d/mu + theta0 <- theta + sumsq0 <- sumsq + e <- l$fn; grad <- l$grad + break + } else{ # increase damping coefficient by a factor of mu + d <- d*mu + } + } + if(termPar < tol) { + WhyStop <- "Tolerance" + break + } + if(abs(sumSqDiff) < 0.01*sumsq0){ + WhyStop <- "No significant change" + break + } + if(i == maxIter){ + WhyStop <- "Maximum Iteration Limit" + break + } + } + # theta <- theta0 + sigma2 <- sum(e^2)/df + vcov <- 1/df*Hinv*sigma2 + + list(params=theta,residuals=e,vcov=vcov,sigma = sqrt(sigma2), + termination=list(WhyStop=WhyStop,iter=i,FcnCount = countObj, + CostFcn=sumsq0)) +} + +#' Create optimization options +#' +#' Specify optimization options that are to be passed to the +#' numerical estimation routines +#' +#' @param tol Minimum 2-norm of the gradient (Default: \code{1e-2}) +#' @param maxIter Maximum number of iterations to be performed +#' @param LMinit Starting value of search-direction length +#' in the Levenberg-Marquardt method (Default: \code{0.01}) +#' @param LMstep Size of the Levenberg-Marquardt step (Default: \code{2}) +#' @param display Argument whether to display iteration details or not +#' (Default: \code{"off"}) +#' +#' @export +optimOptions <- function(tol=1e-2,maxIter=20,LMinit=0.01,LMstep=2, + display=c("off","on")[1]){ + return(list(tol=tol,maxIter= maxIter, + adv= list(LMinit=LMinit,LMstep=LMstep),display=display)) +} + +#' Parameter covariance of the identified model +#' +#' Obtain the parameter covariance matrix of the linear, identified +#' parametric model +#' +#' @param sys a linear, identified parametric model +#' +#' @export +getcov <- function(sys){ + sys$stats$vcov +} + +armaxGrad <- function(theta,e,dots){ + y <- dots[[1]]; u <- dots[[2]]; order <- dots[[3]]; + na <- order[1];nb <- order[2]; nc <- order[3]; nk <- order[4] + nb1 <- nb+nk-1 ; n <- max(na,nb1,nc);N <- dim(y)[1] + + l <- list() + if(is.null(e)){ + e <- dots[[4]]; l$e <- e + } + + yout <- apply(y,2,padZeros,n=n) + uout <- apply(u,2,padZeros,n=n) + eout <- apply(e,2,padZeros,n=n) + + reg <- function(i) { + if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1 + matrix(c(-yout[i-1:na,],uout[v,],eout[i-1:nc,])) + } + + X <- t(sapply(n+1:(N+n),reg)) + Y <- yout[n+1:(N+n),,drop=F] + fn <- Y-X%*%theta + + # Compute Gradient + filt1 <- signal::Arma(b=1,a=c(1,theta[(na+nb+1:nc)])) + grad <- apply(X,2,signal::filter,filt=filt1) + + l$grad <- grad[1:N,,drop=F];l$fn <- fn[1:N,,drop=F] + return(l) +} + +oeGrad <- function(theta,e,dots){ + y <- dots[[1]]; u <- dots[[2]]; order <- dots[[3]]; + nb <- order[1];nf <- order[2]; nk <- order[3]; + nb1 <- nb+nk-1 ; n <- max(nb1,nf) + N <- dim(y)[1] + l <- list() + if(is.null(e)){ + iv <- dots[[4]] + fn <- y-iv; l$e <- fn + } else{ + iv <- y-e + } + uout <- apply(u,2,leftPadZeros,n=n) + ivout <- apply(iv,2,leftPadZeros,n=n) + + reg <- function(i) { + if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1 + matrix(c(uout[v,],-ivout[i-1:nf,])) + } + # Compute new regressor matrix and residuals + X <- t(sapply(n+1:N,reg)) + fn <- y-X%*%theta + + # Compute gradient + filt1 <- signal::Arma(b=1,a=c(1,theta[nb+1:nf,])) + grad <- apply(X,2,signal::filter,filt=filt1) + + l$fn <- fn; l$grad<-grad + return(l) +} + +bjGrad <- function(theta,e,dots){ + y <- dots[[1]]; u <- dots[[2]]; order <- dots[[3]]; + nb <- order[1];nc <- order[2]; nd <- order[3]; + nf <- order[4]; nk <- order[5];nb1 <- nb+nk-1 ; n <- max(nb1,nc,nd,nf); + N <- dim(y)[1] + + l <- list() + if(is.null(e)){ + zeta <- dots[[4]] + w <- y-zeta + e <- dots[[5]]; l$e <- e + } else{ + filt_ts <- signal::Arma(b=c(1,theta[nb+1:nc]), + a=c(1,theta[nb+nc+1:nd])) + w <- matrix(signal::filter(filt_ts,e)) + zeta <- y-w + } + + uout <- apply(u,2,leftPadZeros,n=n) + zetaout <- apply(zeta,2,leftPadZeros,n=n) + eout <- apply(e,2,leftPadZeros,n=n) + wout <- apply(w,2,leftPadZeros,n=n) + + reg <- function(i) { + if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1 + ereg <- if(nc==0) NULL else eout[i-1:nc,] + matrix(c(uout[v,],ereg,wout[i-1:nd,],-zetaout[i-1:nf,])) + } + + # Compute new regressor matrix and residuals + X <- t(sapply(n+1:N,reg)) + fn <- y-X%*%theta + + # Computing gradient + C_params <- if(nc==0) NULL else theta[nb+1:nc] + den <- as.numeric(polynom::polynomial(c(1,C_params))* + polynom::polynomial(c(1,theta[nb+nc+nd+1:nf]))) + filt1 <- signal::Arma(b=c(1,theta[nb+nc+1:nd]), + a=den) + grad <- apply(X,2,signal::filter,filt=filt1) + + l$fn <- fn; l$grad <- grad + return(l) +} + +checkInitSys <- function(init_sys){ + z <- strsplit(toString(sys.call(which=-1)),split = ",")[[1]][1] + if(init_sys$type!=z){ + errMes <- paste("An idpoly model of",toupper(z),"structure expected for the",z,"command.") + stop(errMes) + } +} + +leftPadZeros <- function(x,n) c(rep(0,n),x) +padZeros <- function(x,n) c(rep(0,n),x,rep(0,n)) + +integfilter <- function(x){ + as.numeric(stats::filter(x,filter=c(1,-1),"convolution",sides = 1, + circular = T)) +}
\ No newline at end of file diff --git a/R/estpoly.R b/R/estpoly.R new file mode 100644 index 0000000..e37ad1d --- /dev/null +++ b/R/estpoly.R @@ -0,0 +1,636 @@ +#' Estimated polynomial object +#' +#' Estimated discrete-time polynomial model returned from an estimation +#' routine. +#' +#' @param sys an \code{idpoly} object containing the estimated polynomial +#' coefficients +#' @param fitted.values 1-step ahead predictions on the training dataset +#' @param residuals 1-step ahead prediction errors +#' @param options optimization specification ser used (applicable for non-linear least +#' squares) +#' @param call the matched call +#' @param stats a list containing estimation statistics +#' @param termination termination criteria for optimization +#' @param input input signal of the training data-set +#' +#' @details +#' Do not use \code{estpoly} for directly specifing an input-output polynomial model. +#' \code{\link{idpoly}} is to be used instead +#' +#' @export +estpoly <- function(sys,fitted.values,residuals,options=NULL, + call,stats,termination=NULL,input){ + out <- list(sys=sys,fitted.values=fitted.values, + residuals=residuals,input=input,call=call, + stats=stats,options=options,termination=termination) + class(out) <- "estpoly" + out +} + + +#' @export +print.estpoly <- function(x,...){ + print(summary(x),...) +} + +#' @export +summary.estpoly <- function(object,...) +{ + model <- object$sys + coefs <- params(model) + + se <- sqrt(diag(getcov(object))) + params <- data.frame(Estimated=coefs,se=se) + + report <- list(fit=fitch(object),params=params) + res <- list(model=model,report=report) + class(res) <- "summary.estpoly" + res +} + +#' Fit Characteristics +#' +#' Returns quantitative assessment of the estimated model as a list +#' +#' @param x the estimated model +#' +#' @return +#' A list containing the following elements +#' +#' \item{MSE}{Mean Square Error measure of how well the response of the model fits +#' the estimation data} +#' \item{FPE}{Final Prediction Error} +#' \item{FitPer}{Normalized root mean squared error (NRMSE) measure of how well the +#' response of the model fits the estimation data, expressed as a percentage.} +#' \item{AIC}{Raw Akaike Information Citeria (AIC) measure of model quality} +#' \item{AICc}{Small sample-size corrected AIC} +#' \item{nAIC}{Normalized AIC} +#' \item{BIC}{Bayesian Information Criteria (BIC)} +#' +#' @export +fitch <- function(x){ + y <- fitted(x) + resid(x) + ek <- as.matrix(resid(x)) + N <- nrow(ek); np <- length(params(x$sys)) + + # fit characteristics + mse <- det(t(ek)%*%ek)/N + fpe <- mse*(1+np/N)/(1-np/N) + nrmse <- 1 - sqrt(sum(ek^2))/sqrt(sum((y-mean(y))^2)) + AIC <- N*log(mse) + 2*np + N*dim(matrix(y))[2]*(log(2*pi)+1) + AICc <- AIC*2*np*(np+1)/(N-np-1) + nAIC <- log(mse) + 2*np/N + BIC <- N*log(mse) + N*dim(matrix(y))[2]*(log(2*pi)+1) + np*log(N) + + list(MSE=mse,FPE=fpe,FitPer = nrmse*100,AIC=AIC,AICc=AICc,nAIC=nAIC,BIC=BIC) +} + +#' @export +print.summary.estpoly <- function(x,digits=4,...){ + print(x$model,se=x$report$params[,2],dig=digits) + cat("\n Fit Characteristics \n") + print(data.frame(x$report$fit),digits=digits) +} + +#' @import ggplot2 +#' @export +plot.estpoly <- function(x,newdata=NULL,...){ + + if(is.null(newdata)){ + ypred <- ts(fitted(x),names="Predicted") + yact <- ts(fitted(x) + resid(x),names="Actual") + time <- time(x$input) + titstr <- "Predictions of Model on Training Set" + } else{ + if(class(newdata)!="idframe") stop("Only idframe objects allowed") + ypred <- predict(x,newdata) + yact <- outputData(newdata)[,1] + time <- time(newdata) + titstr <- "Predictions of Model on Test Set" + } + df <- data.frame(Predicted=ypred,Actual=yact,Time=time) + with(df,ggplot(df, aes(x = Actual,y=Predicted)) + ggtitle(titstr) + + geom_abline(intercept=0,slope=1,colour="#D55E00") + geom_point()) +} + +#' Plot residual characteristics +#' +#' Computes the 1-step ahead prediction errors (residuals) for an estimated polynomial +#' model, and plots auto-correlation of the residuals and the +#' cross-correlation of the residuals with the input signals. +#' +#' @param model estimated polynomial model +#' @param newdata an optional dataset on which predictions are to be computed. If +#' not supplied, predictions are computed on the training dataset. +#' @export +residplot <- function(model,newdata=NULL){ + if(is.null(newdata)){ + e <- resid(model); u <- model$input + } else{ + if(class(newdata)!="idframe") stop("Only idframe objects allowed") + e <- newdata$output[,1] - predict(model,newdata)[,1] + u <- newdata$input + } + e <- matrix(e) + acorr <- acf(e[,],plot = F); ccorr <- ccf(u[,1],e[,],plot = F) + par(mfrow=c(2,1),mar=c(3,4,3,2)) + plot(acorr,ci=0.99,main="ACF of residuals") + plot(ccorr,ci=0.99,main="CCF between the input and residuals",ylab="CCF") +} + +#' Estimate ARX Models +#' +#' Fit an ARX model of the specified order given the input-output data +#' +#' @param x an object of class \code{idframe} +#' @param order Specification of the orders: the three integer components +#' (na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) and +#' the input-output delay +#' @param lambda Regularization parameter(Default=\code{0.1}) +#' @param intNoise Logical variable indicating whether to add integrators in +#' the noise channel (Default=\code{FALSE}) +#' @param fixed list containing fixed parameters. If supplied, only \code{NA} entries +#' will be varied. Specified as a list of two vectors, each containing the parameters +#' of polynomials A and B respectively. +#' +#' @details +#' SISO ARX models are of the form +#' \deqn{ +#' y[k] + a_1 y[k-1] + \ldots + a_{na} y[k-na] = b_{nk} u[k-nk] + +#' \ldots + b_{nk+nb} u[k-nk-nb] + e[k] +#' } +#' The function estimates the coefficients using linear least squares (with +#' regularization). +#' \cr +#' The data is expected to have no offsets or trends. They can be removed +#' using the \code{\link{detrend}} function. +#' \cr +#' To estimate finite impulse response(\code{FIR}) models, specify the first +#' order to be zero. +#' +#' @return +#' An object of class \code{estpoly} containing the following elements: +#' \item{sys}{an \code{idpoly} object containing the +#' fitted ARX coefficients} +#' \item{fitted.values}{the predicted response} +#' \item{residuals}{the residuals} +#' \item{input}{the input data used} +#' \item{call}{the matched call} +#' \item{stats}{A list containing the following fields: \cr +#' \code{vcov} - the covariance matrix of the fitted coefficients \cr +#' \code{sigma} - the standard deviation of the innovations\cr +#' \code{df} - the residual degrees of freedom} +#' +#' +#' @references +#' Arun K. Tangirala (2015), \emph{Principles of System Identification: +#' Theory and Practice}, CRC Press, Boca Raton. Section 21.6.1 +#' +#' Lennart Ljung (1999), \emph{System Identification: Theory for the User}, +#' 2nd Edition, Prentice Hall, New York. Section 10.1 +#' +#' @examples +#' data(arxsim) +#' mod_arx <- arx(arxsim,c(1,2,2)) +#' mod_arx +#' plot(mod_arx) # plot the predicted and actual responses +#' +#' @export +arx <- function(x,order=c(1,1,1),lambda=0.1,intNoise=FALSE, + fixed=NULL){ + y <- outputData(x); u <- inputData(x) + if(intNoise){ + y <- apply(y,2,diff) + u <- apply(u,2,diff) + } + N <- dim(y)[1] + na <- order[1];nb <- order[2]; nk <- order[3] + nb1 <- nb+nk-1 ; n <- max(na,nb1); df <- N-na-nb + + yout <- apply(y,2,padZeros,n=n); + uout <- apply(u,2,padZeros,n=n); + + fixedflag = is.null(fixed) + uindex <- nk:nb1 + if(na!=0) yindex <- 1:na + + if(!fixedflag){ + # checking for correct specification of fixed parameters + fixedA <- NULL;fixedB <- NULL + g(fixedA,fixedB) %=% lapply(fixed,length) + if(fixedA != na || fixedB != nb) + stop("Number of parameters incorrectly specified in 'fixed'") + + fixedpars <- unlist(fixed) + fixedpars <- fixedpars[!is.na(fixedpars)] + df <- df + length(fixedpars) + + fixedpos_B <- which(!is.na(fixed[[2]])) + uindex <- uindex[!uindex %in% (nk+fixedpos_B-1)] + fixedpos_A <- numeric(0) + if(na!=0){ + fixedpos_A <- which(!is.na(fixed[[1]])) + yindex <- yindex[!yindex %in% fixedpos_A] + } + } + + reg <- function(i) { + # regressor + temp <- numeric(0) + if(na!=0) temp <- c(temp,-yout[i-yindex,]) + phi <- t(c(temp,uout[i-uindex,])) + l <- list(phi=phi) + # fixed regressors + if(!fixedflag){ + temp <- numeric(0) + if(length(fixedpos_A)!=0){ + temp <- c(temp,-yout[i-fixedpos_A,]) + } + if(length(fixedpos_B)!=0){ + temp <- c(temp,uout[i-(nk+fixedpos_B),]) + } + l$fixed <- t(temp) + } + return(l) + } + + temp <- lapply(n+1:(N+n),reg) + X <- do.call(rbind,lapply(temp, function(x) x[[1]])) + Y <- yout[n+1:(N+n),,drop=F] + fixedY <- matrix(rep(0,nrow(Y))) + if(!fixedflag){ + fixedreg <- do.call(rbind,lapply(temp, function(x) x[[2]])) + fixedY <- fixedreg%*%fixedpars + Y <- Y - fixedY + } + + # lambda <- 0.1 + inner <- t(X)%*%X + lambda*diag(dim(X)[2]) + innerinv <- solve(inner) + pinv <- innerinv%*% t(X) + coef <- pinv%*%Y + + eps <- X%*%coef + sigma2 <- sum(eps^2)/(df+n) + vcov <- sigma2 * innerinv + fit <- (X%*%coef+fixedY)[1:N,,drop=F] + if(intNoise) fit <- apply(fit,2,cumsum) + + if(!fixedflag){ + findex <- which(!is.na(unlist(fixed))) + eindex <- 1:(na+nb) + eindex <- eindex[!eindex %in% findex] + temp <- rep(0,na+nb); temp2 <- matrix(0,nrow=na+nb,ncol=na+nb) + + temp[eindex] <- coef; temp[findex] <- fixedpars; coef <- temp + temp2[eindex,eindex] <- vcov; vcov <- temp2 + } + + if(na==0){ + A <- 1 + } else { + A <- c(1,coef[1:na]) + } + model <- idpoly(A = A,B = coef[na+1:nb], + ioDelay = nk,Ts=deltat(x),noiseVar = sqrt(sigma2), + intNoise=intNoise,unit=x$unit) + + estpoly(sys = model,stats=list(vcov = vcov, sigma = sqrt(sigma2), + df = df),fitted.values=fit,residuals=eps[1:N,,drop=F], + call=match.call(),input=u) +} + +#' Estimate ARMAX Models +#' +#' Fit an ARMAX model of the specified order given the input-output data +#' +#' @param x an object of class \code{idframe} +#' @param order Specification of the orders: the four integer components +#' (na,nb,nc,nk) are the order of polynolnomial A, order of polynomial B +#' + 1, order of the polynomial C,and the input-output delay respectively +#' @param init_sys Linear polynomial model that configures the initial parameterization. +#' Must be an ARMAX model. Overrules the \code{order} argument +#' @param intNoise Logical variable indicating whether to add integrators in +#' the noise channel (Default=\code{FALSE}) +#' @param options Estimation Options, setup using \code{\link{optimOptions}} +#' +#' @details +#' SISO ARMAX models are of the form +#' \deqn{ +#' y[k] + a_1 y[k-1] + \ldots + a_{na} y[k-na] = b_{nk} u[k-nk] + +#' \ldots + b_{nk+nb} u[k-nk-nb] + c_{1} e[k-1] + \ldots c_{nc} e[k-nc] +#' + e[k] +#' } +#' The function estimates the coefficients using non-linear least squares +#' (Levenberg-Marquardt Algorithm) +#' \cr +#' The data is expected to have no offsets or trends. They can be removed +#' using the \code{\link{detrend}} function. +#' +#' @return +#' An object of class \code{estpoly} containing the following elements: +#' \item{sys}{an \code{idpoly} object containing the +#' fitted ARMAX coefficients} +#' \item{fitted.values}{the predicted response} +#' \item{residuals}{the residuals} +#' \item{input}{the input data used} +#' \item{call}{the matched call} +#' \item{stats}{A list containing the following fields: \cr +#' \code{vcov} - the covariance matrix of the fitted coefficients \cr +#' \code{sigma} - the standard deviation of the innovations} +#' \item{options}{Option set used for estimation. If no +#' custom options were configured, this is a set of default options} +#' \item{termination}{Termination conditions for the iterative +#' search used for prediction error minimization: +#' \code{WhyStop} - Reason for termination \cr +#' \code{iter} - Number of Iterations \cr +#' \code{iter} - Number of Function Evaluations } +#' +#' +#' @references +#' Arun K. Tangirala (2015), \emph{Principles of System Identification: +#' Theory and Practice}, CRC Press, Boca Raton. Sections 14.4.1, 21.6.2 +#' +#' @examples +#' data(armaxsim) +#' z <- dataSlice(armaxsim,end=1533) # training set +#' mod_armax <- armax(z,c(1,2,1,2)) +#' mod_armax +#' +#' @export +armax <- function(x,order=c(0,1,1,0),init_sys=NULL,intNoise=FALSE, + options=optimOptions()){ + y <- outputData(x); u <- inputData(x) + if(intNoise){ + y <- apply(y,2,diff) + u <- apply(u,2,diff) + } + N <- dim(y)[1] + if(!is.null(init_sys)){ + checkInitSys(init_sys) + + # Extract orders from initial guess + na <- length(init_sys$A) -1;nb <- length(init_sys$B); + nc <- length(init_sys$C) -1;nk <- init_sys$ioDelay + order <- c(na,nb,nc,nk) + + # Initial guess + theta0 <- matrix(params(init_sys)) + ivs <- matrix(predict(init_sys,x)) + e_init <- y-ivs + } else{ + na <- order[1];nb <- order[2]; nc <- order[3]; nk <- order[4] + + if(nc<1) + stop("Error: Not an ARMAX model") + + # Initial Parameter Estimates + mod_arx <- iv4(x,c(na,nb,nk)) # fitting ARX model + eps_init <- matrix(resid(mod_arx)) + mod_ma <- arima(eps_init,order=c(0,0,nc),include.mean = F) + e_init <- matrix(mod_ma$residuals); e_init[is.na(e_init)] <- 0 + theta0 <- matrix(c(mod_arx$sys$A[-1],mod_arx$sys$B,mod_ma$coef)) + } + + nb1 <- nb+nk-1 ; n <- max(na,nb1,nc); df <- N - na - nb - nc + l <- levbmqdt(y,u,order,e_init,obj=armaxGrad, + theta0=theta0,N=N,opt=options) + theta <- l$params + e <- ts(l$residuals,start = start(y),deltat = deltat(y)) + fit <- matrix(y-e) + if(intNoise) fit <- apply(fit,2,cumsum) + + model <- idpoly(A = c(1,theta[1:na]),B = theta[na+1:nb], + C = c(1,theta[na+nb+1:nc]),ioDelay = nk,Ts=deltat(x), + noiseVar = l$sigma,intNoise=intNoise,unit=x$unit) + + estpoly(sys = model,stats=list(vcov = l$vcov, sigma = l$sigma), + fitted.values=fit,residuals=e,call=match.call(),input=u, + options = options,termination = l$termination) +} + +#' Estimate Output-Error Models +#' +#' Fit an output-error model of the specified order given the input-output data +#' +#' @param x an object of class \code{idframe} +#' @param order Specification of the orders: the four integer components +#' (nb,nf,nk) are order of polynomial B + 1, order of the polynomial F, +#' and the input-output delay respectively +#' @param init_sys Linear polynomial model that configures the initial parameterization. +#' Must be an OE model. Overrules the \code{order} argument +#' @param options Estimation Options, setup using +#' \code{\link{optimOptions}} +#' +#' @details +#' SISO OE models are of the form +#' \deqn{ +#' y[k] + f_1 y[k-1] + \ldots + f_{nf} y[k-nf] = b_{nk} u[k-nk] + +#' \ldots + b_{nk+nb} u[k-nk-nb] + f_{1} e[k-1] + \ldots f_{nf} e[k-nf] +#' + e[k] +#' } +#' The function estimates the coefficients using non-linear least squares +#' (Levenberg-Marquardt Algorithm) +#' \cr +#' The data is expected to have no offsets or trends. They can be removed +#' using the \code{\link{detrend}} function. +#' +#' @return +#' An object of class \code{estpoly} containing the following elements: +#' \item{sys}{an \code{idpoly} object containing the +#' fitted OE coefficients} +#' \item{fitted.values}{the predicted response} +#' \item{residuals}{the residuals} +#' \item{input}{the input data used} +#' \item{call}{the matched call} +#' \item{stats}{A list containing the following fields: \cr +#' \code{vcov} - the covariance matrix of the fitted coefficients \cr +#' \code{sigma} - the standard deviation of the innovations} +#' \item{options}{Option set used for estimation. If no +#' custom options were configured, this is a set of default options} +#' \item{termination}{Termination conditions for the iterative +#' search used for prediction error minimization: +#' \code{WhyStop} - Reason for termination \cr +#' \code{iter} - Number of Iterations \cr +#' \code{iter} - Number of Function Evaluations } +#' +#' @references +#' Arun K. Tangirala (2015), \emph{Principles of System Identification: +#' Theory and Practice}, CRC Press, Boca Raton. Sections 14.4.1, 17.5.2, +#' 21.6.3 +#' +#' @examples +#' data(oesim) +#' z <- dataSlice(oesim,end=1533) # training set +#' mod_oe <- oe(z,c(2,1,2)) +#' mod_oe +#' plot(mod_oe) # plot the predicted and actual responses +#' +#' @export +oe <- function(x,order=c(1,1,0),init_sys=NULL,options=optimOptions()){ + y <- outputData(x); u <- inputData(x); N <- dim(y)[1] + if(!is.null(init_sys)){ + checkInitSys(init_sys) + + # Extract orders from initial guess + nb <- length(init_sys$B); nf <- length(init_sys$F1) -1 + nk <- init_sys$ioDelay;order <- c(nb,nf,nk) + + # Initial guess + theta0 <- matrix(params(init_sys)) + ivs <- matrix(predict(init_sys,x)) + e_init <- y-ivs + } else{ + nb <- order[1];nf <- order[2]; nk <- order[3]; + nb1 <- nb+nk-1 ; n <- max(nb1,nf); + + if(nf<1) + stop("Not an OE model") + + # Initial Model + mod_arx <- iv4(x,c(nf,nb,nk)) # fitting ARX model + wk <- resid(mod_arx) + e_init <- as.numeric(stats::filter(wk,filter=-mod_arx$sys$A[-1], + method = "recursive")) + ivs <- y-e_init + theta0 <- matrix(c(mod_arx$sys$B,mod_arx$sys$A[-1])) + } + nb1 <- nb+nk-1 ; n <- max(nb1,nf);df <- N - nb - nf + + l <- levbmqdt(y,u,order,ivs,obj=oeGrad,theta0=theta0,N=N, + opt=options) + theta <- l$params + e <- ts(l$residuals,start = start(y),deltat = deltat(y)) + + model <- idpoly(B = theta[1:nb],F1 = c(1,theta[nb+1:nf]), + ioDelay = nk,Ts=deltat(x),noiseVar = l$sigma,unit=x$unit) + + estpoly(sys = model,stats=list(vcov = l$vcov, sigma = l$sigma), + fitted.values=y-e,residuals=e,call=match.call(),input=u, + options = options,termination = l$termination) +} + +#' Estimate Box-Jenkins Models +#' +#' Fit a box-jenkins model of the specified order from input-output data +#' +#' @param z an \code{idframe} object containing the data +#' @param order Specification of the orders: the five integer components +#' (nb,nc,nd,nf,nk) are order of polynomial B + 1, order of the polynomial C, +#' order of the polynomial D, order of the polynomial F, and the +#' input-output delay respectively +#' @param init_sys Linear polynomial model that configures the initial parameterization. +#' Must be a BJ model. Overrules the \code{order} argument +#' @param options Estimation Options, setup using +#' \code{\link{optimOptions}} +#' +#' @details +#' SISO BJ models are of the form +#' \deqn{ +#' y[k] = \frac{B(q^{-1})}{F(q^{-1})}u[k-nk] + +#' \frac{C(q^{-1})}{D(q^{-1})} e[k] +#' } +#' The orders of Box-Jenkins model are defined as follows: +#' \deqn{ +#' B(q^{-1}) = b_1 + b_2q^{-1} + \ldots + b_{nb} q^{-nb+1} +#' } +#' +#' \deqn{ +#' C(q^{-1}) = 1 + c_1q^{-1} + \ldots + c_{nc} q^{-nc} +#' } +#' +#' \deqn{ +#' D(q^{-1}) = 1 + d_1q^{-1} + \ldots + d_{nd} q^{-nd} +#' } +#' \deqn{ +#' F(q^{-1}) = 1 + f_1q^{-1} + \ldots + f_{nf} q^{-nf} +#' } +#' +#' The function estimates the coefficients using non-linear least squares +#' (Levenberg-Marquardt Algorithm) +#' \cr +#' The data is expected to have no offsets or trends. They can be removed +#' using the \code{\link{detrend}} function. +#' +#' @return +#' An object of class \code{estpoly} containing the following elements: +#' \item{sys}{an \code{idpoly} object containing the +#' fitted BJ coefficients} +#' \item{fitted.values}{the predicted response} +#' \item{residuals}{the residuals} +#' \item{input}{the input data used} +#' \item{call}{the matched call} +#' \item{stats}{A list containing the following fields: \cr +#' \code{vcov} - the covariance matrix of the fitted coefficients \cr +#' \code{sigma} - the standard deviation of the innovations} +#' \item{options}{Option set used for estimation. If no +#' custom options were configured, this is a set of default options} +#' \item{termination}{Termination conditions for the iterative +#' search used for prediction error minimization: +#' \code{WhyStop} - Reason for termination \cr +#' \code{iter} - Number of Iterations \cr +#' \code{iter} - Number of Function Evaluations } +#' +#' @references +#' Arun K. Tangirala (2015), \emph{Principles of System Identification: +#' Theory and Practice}, CRC Press, Boca Raton. Sections 14.4.1, 17.5.2, +#' 21.6.3 +#' +#' @examples +#' data(bjsim) +#' z <- dataSlice(bjsim,end=1500) # training set +#' mod_bj <- bj(z,c(2,1,1,1,2)) +#' mod_bj +#' residplot(mod_bj) # residual plots +#' +#' @export +bj <- function(z,order=c(1,1,1,1,0), + init_sys=NULL,options=optimOptions()){ + y <- outputData(z); u <- inputData(z); N <- dim(y)[1] + if(!is.null(init_sys)){ + checkInitSys(init_sys) + + # Extract orders from initial guess + nb <- length(init_sys$B); nf <- length(init_sys$F1) -1 + nc <- length(init_sys$C) -1;nd <- length(init_sys$d) -1 + nk <- init_sys$ioDelay;order <- c(nb,nc,nd,nf,nk) + + # Initial guess + theta0 <- matrix(params(init_sys)) + ivs <- matrix(predict(init_sys,z)) + e_init <- y-ivs + } else{ + nb <- order[1];nc <- order[2]; nd <- order[3]; + nf <- order[4]; nk <- order[5]; + + if(nc==0 && nd==0){ + oe(z,c(nb,nf,nk)) + } else{ + + # Initial Guess + mod_oe <- oe(z,c(nb,nf,nk)) + v <- resid(mod_oe); zeta <- matrix(predict(mod_oe)) + mod_arma <- arima(v,order=c(nd,0,nc),include.mean = F) + C_params <- if(nc==0) NULL else coef(mod_arma)[nd+1:nc] + theta0 <- matrix(c(mod_oe$sys$B,C_params, + -coef(mod_arma)[1:nd],mod_oe$sys$F1[-1])) + eps <- matrix(resid(mod_arma)) + } + } + l <- levbmqdt(y,u,order,zeta,eps,obj=bjGrad,theta0=theta0,N=N, + opt=options) + theta <- l$params + e <- ts(l$residuals,start = start(y),deltat = deltat(y)) + + C_params <- if(nc==0) NULL else theta[nb+1:nc] + model <- idpoly(B = theta[1:nb],C=c(1,C_params), + D=c(1,theta[nb+nc+1:nd]), + F1 = c(1,theta[nb+nc+nd+1:nf]), + ioDelay = nk,Ts=deltat(z),noiseVar = l$sigma,unit=z$unit) + + estpoly(sys = model,stats=list(vcov = l$vcov, sigma = l$sigma), + fitted.values=y-e,residuals=e,call=match.call(),input=u, + options = options,termination = l$termination) + +}
\ No newline at end of file diff --git a/R/idframe.R b/R/idframe.R new file mode 100644 index 0000000..e619fca --- /dev/null +++ b/R/idframe.R @@ -0,0 +1,280 @@ +#' S3 class for storing input-output data. +#' +#' \code{idframe} is an S3 class for storing and manipulating input-ouput data. It supports discrete time and frequency domain data. +#' +#' @param output dataframe/matrix/vector containing the outputs +#' @param input dataframe/matrix/vector containing the inputs +#' @param Ts sampling interval (Default: 1) +#' @param start Time of the first observation +#' @param end Time of the last observation Optional Argument +#' @param unit Time unit (Default: "seconds") +#' +#' @return an idframe object +#' +#' @seealso \code{\link{plot.idframe}}, the plot method for idframe objects +#' +#' @examples +#' +#' dataMatrix <- matrix(rnorm(1000),ncol=5) +#' data <- idframe(output=dataMatrix[,3:5],input=dataMatrix[,1:2],Ts=1) +#' +#' @export +idframe <- function(output,input=NULL,Ts = 1,start=0,end=NULL, + unit = c("seconds","minutes","hours", + "days")[1]){ + + l <- list(output) + if(!is.null(input)) l[[2]] <- input + l2 <- lapply(l,data.frame) + n <- dim(l2[[1]]) + dims <- sapply(l2,ncol) + colnames(l2[[1]]) <- sapply(1:dims[1], + function(x) paste("y",x,sep = "")) + if(!is.null(input)){ + colnames(l2[[2]]) <- sapply(1:dims[2], + function(x) paste("u",x,sep = "")) + } + + if(!is.null(end)){ + start <- end - Ts*(n-1) + } + + l3 <- lapply(l2,ts,start=start,deltat=Ts) + + if(!is.null(input)) input <- l3[[2]] + # Object Constructor + dat <- list(output=l3[[1]],input=input,unit=unit) + class(dat) <- "idframe" + return(dat) +} + +#' Plotting idframe objects +#' +#' Plotting method for objects inherting from class \code{idframe} +#' +#' @param x an \code{idframe} object +#' @param col line color, to be passed to plot.(Default=\code{"steelblue"}) +#' @param lwd line width, in millimeters(Default=\code{1}) +#' @param main the plot title. (Default = \code{NULL}) +#' @param size text size (Default = \code{12}) +#' @param \ldots additional arguments +#' +#' @examples +#' data(cstr) +#' plot(cstr,col="blue") +#' +#' @import ggplot2 reshape2 +#' +#' @export +plot.idframe <- function(x,col="steelblue",lwd=1,main=NULL,size=12,...){ + loadNamespace("ggplot2") + if(nInputSeries(x)==0){ + data <- outputData(x) + } else{ + data <- cbind(outputData(x),inputData(x)) + colnames(data) <- c(outputNames(x),inputNames(x)) + } + df <- melt(data.frame(time=as.numeric(time(data)), data), id.vars="time") + with(df,ggplot(df, aes(time, value)) + geom_line(size=lwd,color=col) + + facet_grid(variable ~ .,scales="free") + theme_bw(size) + + ylab("Amplitude") + ggtitle(main) + + xlab(paste("Time (",x$unit,")",sep=""))) + #theme(axis.title.x=element_text(size=11)) + ggtitle(main) +} + +#' @export +summary.idframe <- function(object,...){ + out_sum <- summary(outputData(object)) + in_sum <- summary(inputData(object)) + + out <- list(out_sum=out_sum,in_sum=in_sum,Ts=deltat(object), + unit=object$unit,nsample = dim(outputData(object))[1]) + + class(out) <- "summary.idframe" + return(out) +} + +#' @export +print.summary.idframe <- function(x,...){ + cat("\t\t Number of samples:");cat(x$nsample) + cat("\nSampling time: ") + cat(x$Ts);cat(" ");cat(x$unit) + + cat("\n\n") + cat("Outputs \n") + print(x$out_sum) + cat("\n") + + cat("Inputs \n") + print(x$in_sum) +} + +#' Sampling times of IO data +#' +#' \code{time} creates the vector of times at which data was sampled. \code{frequency} returns the number of damples per unit time and \code{deltat} the time-interval +#' between observations +#' +#' @param x a idframe object, or a univariate or multivariate time-series, or a vector or matrix +#' +#' @aliases frequency deltat +#' +#' @export +time <- function(x){ + if(class(x)[1]!="idframe"){ + stats::time(x) + } else{ + stats::time(x$output) + } +} + +#' @export +frequency <- function(x){ + if(class(x)[1]!="idframe"){ + stats::frequency(x) + } else{ + stats::frequency(x$output) + } +} + +#' @export +deltat <- function(x){ + if(class(x)[1]!="idframe"){ + stats::deltat(x) + } else{ + stats::deltat(x$output) + } +} + +#' S3 class constructor for storing frequency response data +#' +#' @param respData frequency response data. For SISO systems, supply a +#' vector of frequency response values. For MIMO systems with Ny +#' outputs and Nu inputs, supply an array of size c(Ny,Nu,Nw). +#' @param freq frequency points of the response +#' @param Ts sampling time of data +#' @param spec power spectra and cross spectra of the system +#' output disturbances (noise). Supply an array of size (Ny,Ny,Nw) +#' @param covData response data covariance matrices. Supply an array +#' of size (Ny,Nu,Nw,2,2). covData[ky,ku,kw,,] is the covariance matrix +#' of respData[ky,ku,kw] +#' @param noiseCov power spectra variance. Supply an array of +#' size (Ny,Ny,Nw) +#' +#' @return an idfrd object +#' +#' @seealso +#' \code{\link{plot.idfrd}} for generating bode plots, +#' \code{\link{spa}} and \code{\link{etfe}} for estimating the +#' frequency response given input/output data +#' +#' @export +idfrd <- function(respData,freq,Ts,spec=NULL,covData=NULL, + noiseCov=NULL){ + # For SISO systems + if(is.vector(respData)){ + dim(respData) <- c(1,1,nrow(freq)) + if(!is.null(spec)) dim(spec) <- c(1,1,nrow(freq)) + } + + out <- list(response=respData,freq=freq,Ts=Ts,spec=spec,covData= + covData,noiseCov = noiseCov) + class(out) <- "idfrd" + return(out) +} + +#' Plotting idfrd objects +#' +#' Generates the bode plot of the given frequency response data. It uses the +#' ggplot2 plotting engine +#' +#' @param x An object of class \code{idframe} +#' @param col a specification for the line colour (Default : \code{" +#' steelblue"}) +#' @param lwd the line width, a positive number, defaulting to 1 +#' @param \ldots additional arguments +#' +#' @seealso \code{\link[ggplot2]{ggplot}} +#' +#' @examples +#' data(frd) +#' plot(frd) +#' +#' @import ggplot2 reshape2 signal +#' +#' @export +plot.idfrd <- function(x,col="steelblue",lwd=1,...){ + loadNamespace("ggplot2") + nfreq <- dim(x$freq)[1] + mag <- 20*log10(Mod(x$resp)) + nout <- dim(mag)[1]; nin <- dim(mag)[2] + dim(mag) <- c(nin*nout,nfreq) + + temp <- aperm(Arg(x$resp),c(3,2,1));dim(temp) <- c(nfreq,nin*nout) + l <- t(split(temp, rep(1:ncol(temp), each = nrow(temp)))) + phase <- 180/pi*t(sapply(l,signal::unwrap)) + + g <- vector("list",nin*nout) + + for(i in 1:length(g)){ + df <- data.frame(Frequency=x$freq[,],Magnitude=mag[i,], + Phase = phase[i,]) + melt_df <- reshape2::melt(df,id.var="Frequency") + yindex <- (i-1)%/%nin + 1;uindex <- i-nin*(yindex-1) + subtitle <- paste("From: u",uindex," to y",yindex,sep="") + g[[i]] <- with(melt_df,ggplot(melt_df, aes(Frequency, value)) + + geom_line(size=lwd,color=col) + scale_x_log10(expand = c(0.01,0.01)) + + facet_grid(variable ~ .,scales="free_y") + + theme_bw(14) + ylab("") + ggtitle(subtitle) + + xlab(ifelse(yindex==nout,"Frequency","")) + + theme(axis.title.x=element_text(color = "black",face = "plain"), + title=element_text(size=9,color = "black",face="bold")) + + geom_vline(xintercept=max(x$freq),size=1)) + } + + multiplot(plotlist=g,layout=matrix(1:length(g),nrow=nout,byrow=T)) +} + +# Multiple plot function +# +# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects) +# - cols: Number of columns in layout +# - layout: A matrix specifying the layout. If present, 'cols' is ignored. +# +# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE), +# then plot 1 will go in the upper left, 2 will go in the upper right, and +# 3 will go all the way across the bottom. +# +multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { + + # Make a list from the ... arguments and plotlist + plots <- c(list(...), plotlist) + + numPlots = length(plots) + + # If layout is NULL, then use 'cols' to determine layout + if (is.null(layout)) { + # Make the panel + # ncol: Number of columns of plots + # nrow: Number of rows needed, calculated from # of cols + layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), + ncol = cols, nrow = ceiling(numPlots/cols)) + } + + if (numPlots==1) { + print(plots[[1]]) + + } else { + # Set up the page + grid::grid.newpage() + grid::pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) + + # Make each plot, in the correct location + for (i in 1:numPlots) { + # Get the i,j matrix positions of the regions that contain this subplot + matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) + + print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, + layout.pos.col = matchidx$col)) + } + } +} diff --git a/R/idinput.R b/R/idinput.R new file mode 100644 index 0000000..46b8cad --- /dev/null +++ b/R/idinput.R @@ -0,0 +1,156 @@ +#' function to generate input singals (rgs/rbs/prbs/sine) +#' +#' \code{idinput} is a function for generating input signals (rgs/rbs/prbs/sine) for identification purposes +#' +#' @param n integer length of the input singal to be generated +#' @param type the type of input signal to be generated. +#' 'rgs' - generates random gaussian signal +#' 'rbs' - generates random binary signal +#' 'prbs' - generates pseudorandom binary signal +#' 'sine' - generates a signal that is a sum of sinusoids +#' +#' Default value is type='rgs' +#' @param band determines the frequency content of the signal. +#' For type='rbs'/'sine'/, band = [wlow,whigh] +#' which specifies the lower and the upper bound of the passband frequencies(expressed as fractions of Nyquist frequency). Default is c(0,1) +#' For type='prbs', band=[0,B] +#' where B is such that the singal is constant over 1/B (clock period). Default is c(0,1) +#' @param levels row vector defining the input level. It is of the form +#' levels=c(minu, maxu) +#' For 'rbs','prbs', 'sine', the generated signal always between minu and maxu. +#' For 'rgs', minu=mean value of signal minus one standard deviation and maxu=mean value of signal plus one standard deviation +#' +#' Default value is levels=c(-1,1) +#' +#' @export +idinput<-function(n,type='rgs',band=c(0,1),levels=c(-1,1)){ + if(type=="rbs"){ + rbs(n,band,levels) + } else if(type=="rgs"){ + rgs(n,band,levels) + } else if(type=="sine"){ + multisine(n,1,band,levels) + }else if(type=="prbs"){ + idin.prbs(n,band,levels) + } +} + +rgs <- function(n,band,levels){ + u <- butter_filt(rnorm(n),band) + mu<-(levels[1]+levels[2])/2 + sigma<-(levels[2]-levels[1])/2 + u*sigma+mu +} + +rbs <- function(n,band,levels){ + u <- butter_filt(rnorm(n),band) + sapply(u,function(x) if(x>0) levels[2] else levels[1]) +} + +butter_filt <- function(x,band){ + filt <- T; type <- "pass" + if(band[1]<=2e-3){ + if(band[2]==1){ + filt <- F + } else{ + type <- "low" + } + } else{ + if(band[2]==1){ + type <- "high" + } + } + if(filt==T){ + if(type=="low"){ + bf <- signal::butter(8,band[2],type,"z") + } else if(type=="pass"){ + bf <- signal::butter(8,band,type,"z") + }else{ + bf <- signal::butter(8,band[1],type,"z") + } + x <- as.numeric(signal::filter(bf,x)) + } + return(matrix(x,ncol=1)) +} + +multisine <- function(N,nin=1,band,levels){ + sinedata <- list(nSin=10,nTrial=10,gridSkip=1) + freq <- 2*pi*seq(1,floor(N/2),by=sinedata$gridSkip)/N + band <- band*pi + freq <- freq[freq>=band[1] & freq<=band[2]] + nl <- length(freq) + freqInd <- seq(from=1,to=nl,length.out = nin*sinedata$nSin) + + # Begin Trials + trials <- function(x){ + freqs <- freq[freqInd[seq(from=x,to=nin*sinedata$nSin,by=nin)]] + for(k in 1:sinedata$nTrial){ + utrial <- rowSums(cos(sapply(freqs, + function(x) (1:N-1)*x+2*pi*rnorm(1)))) + if(k==1) u <- utrial; temp <- diff(range(utrial)) + if(diff(range(utrial))< temp) { + u <- utrial; temp <- diff(range(utrial)) + } + } + clevel <- max(abs(u)) + diff(levels)/2/clevel*(u+clevel)+levels[1] + } + u <- sapply(1:nin,trials) + return(u) +} + +#' @import bitops signal +idin.prbs<-function(n,band=c(0,1),levels=c(0,1)){ + u <- numeric(0) + for(i in 1:18){ + if(n / (2^i) < 1){ + u <- idin.prbs12(i,band,levels) + break + } + } + return(u[1:n]) +} + +idin.prbs12 <- function(N,band=c(0,1),levels=c(0,1)){ + first=ceiling(abs(rnorm(1)*10)) #some non-zero initial state + x= first + v = vector() + n=2^N-1 + i=1 + clock=floor(1/band[2]) + k=1 + M=rbind(c(0,0,0,0),c(1,2,0,0),c(1,3,0,0),c(1,4,0,0),c(2,5,0,0),c(1,6,0,0), + c(1,7,0,0),c(1,2,7,8),c(4,9,0,0),c(3,10,0,0),c(9,11,0,0), + c(6,8,11,12),c(9,10,12,13),c(4,8,13,14),c(14,15,0,0),c(4,13,15,16), + c(14,17,0,0),c(11,18,0,0)) + repeat{ + a=M[N,1] + b=M[N,2] + c=M[N,3] + d=M[N,4] + four=c(8,12,13,14,16) + if(N %in% four){ + e=bitwXor(bitwShiftR(x,N-a),bitwShiftR(x,N-b)) + f=bitwXor(bitwShiftR(x,N-c),bitwShiftR(x,N-d)) + newbit=bitwAnd(bitwXor(e,f),1) + }else{ + newbit=bitwAnd(bitwXor(bitwShiftR(x,N-a),bitwShiftR(x,N-b)),1) + } + if(k>=clock || i==1){ + v[i]=newbit + i=i+1 + k=1 + }else{ + v[i]=v[i-1] + i=i+1 + k=k+1 + } + x=bitwOr(bitwShiftR(x,1),bitwShiftL(newbit,N-1)) + if(x==first){ + break + } + } + + v=sapply(v, function(x){if (x==0) levels[1] else levels[2]}) + return(v) +}
\ No newline at end of file diff --git a/R/imports.R b/R/imports.R new file mode 100644 index 0000000..e994228 --- /dev/null +++ b/R/imports.R @@ -0,0 +1,5 @@ +#' @importFrom stats predict resid arima ts start ar acf ccf coef fitted lm mvfft rnorm window +#' @importFrom graphics par plot +#' @importFrom utils read.table tail +#' @importFrom grid viewport pushViewport grid.layout grid.newpage +NULL
\ No newline at end of file diff --git a/R/ioNamesData.R b/R/ioNamesData.R new file mode 100644 index 0000000..f28f7ae --- /dev/null +++ b/R/ioNamesData.R @@ -0,0 +1,129 @@ +#' Number of series in input or output +#' +#' Number of series in input or output in a idframe object +#' +#' +#' @param data \code{idframe} object +#' +#' @aliases nOutputSeries +#' +#' @export +nInputSeries <- function(data) { + ifelse(is.null(data$input),0,ncol(data$input)) +} + +#' @export +nOutputSeries <- function(data) ncol(data$output) + +#' Output or Input-data +#' +#' Extract output-data or input-data in idframe objects +#' +#' @param x \code{idframe} object +#' @param series the indices to extract +#' +#' @aliases inputData.idframe outputData outputData.idframe +#' +#' @export +inputData <- function(x,series) UseMethod("inputData") + +#' @export +inputData.default <- function(x,series=NULL){ + print("Not defined for this datatype") +} + +#' @import tframe +#' @export +inputData.idframe <- function(x,series=seq(nInputSeries(x))){ + if (is.null(x$input)) + NULL + else selectSeries(x$input,series=series) +} + +"inputData<-" <- function(x,value) UseMethod("inputData<-") + +"inputData<-.idframe" <- function(x,value){ + x$input <- value + x +} + + +#' @export +inputNames <- function(x) UseMethod("inputNames") + +#' @export +inputNames.default <- function(x){ + print("Not defined for this datatype") +} + +#' @import tframe +#' @export +inputNames.idframe <- function(x){ + seriesNames(inputData(x)) +} + +#' Extract or set series' names +#' +#' Extract or set names of series in input or output +#' +#' @param x \code{idframe} object +#' @param value vector of strings +#' +#' @aliases inputNames inputNames<-.idframe outputNames outputNames outputNames<- outputNames<-.idframe +#' +#' @export +"inputNames<-" <- function(x,value) UseMethod("inputNames<-") + +#' @import tframe +#' @export +"inputNames<-.idframe" <- function(x,value){ + seriesNames(inputData(x)) <- value + x +} + +#' @export +outputData <- function(x,series) UseMethod("outputData") + +#' @export +outputData.default <- function(x,series=NULL){ + print("Not defined for this datatype") +} + +#' @import tframe +#' @export +outputData.idframe <- function(x,series=seq(nOutputSeries(x))){ + if (is.null(x$output)) + NULL + else selectSeries(x$output,series=series) +} + +"outputData<-" <- function(x,value) UseMethod("outputData<-") + +"outputData<-.idframe" <- function(x,value){ + x$output <- value + x +} + +#' @export +outputNames <- function(x) UseMethod("outputNames") + +#' @export +outputNames.default <- function(x){ + print("Not defined for this datatype") +} + +#' @import tframe +#' @export +outputNames.idframe <- function(x){ + seriesNames(outputData(x)) +} + +#' @export +"outputNames<-" <- function(x,value) UseMethod("outputNames<-") + +#' @import tframe +#' @export +"outputNames<-.idframe" <- function(x,value){ + seriesNames(outputData(x)) <- value + x +} @@ -0,0 +1,157 @@ +#' ARX model estimation using instrumental variable method +#' +#' Estimates an ARX model of the specified order from input-output data using +#' the instrument variable method. If arbitrary instruments are not supplied +#' by the user, the instruments are generated using the arx routine +#' +#' @param z an idframe object containing the data +#' @param order Specification of the orders: the three integer components +#' (na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) +#' and the input-output delay +#' @param x instrument variable matrix. x must be of the same size as the output +#' data. (Default: \code{NULL}) +#' +#' @return +#' An object of class \code{estpoly} containing the following elements: +#' \item{sys}{an \code{idpoly} object containing the +#' fitted ARX coefficients} +#' \item{fitted.values}{the predicted response} +#' \item{residuals}{the residuals} +#' \item{input}{the input data used} +#' \item{call}{the matched call} +#' \item{stats}{A list containing the following fields: \cr +#' \code{vcov} - the covariance matrix of the fitted coefficients \cr +#' \code{sigma} - the standard deviation of the innovations\cr +#' \code{df} - the residual degrees of freedom} +#' +#' @references +#' Arun K. Tangirala (2015), \emph{Principles of System Identification: +#' Theory and Practice}, CRC Press, Boca Raton. Sections 21.7.1, 21.7.2 +#' +#' Lennart Ljung (1999), \emph{System Identification: Theory for the User}, +#' 2nd Edition, Prentice Hall, New York. Section 7.6 +#' +#' @examples +#' data(arxsim) +#' mod_iv <- iv(arxsim,c(2,1,1)) +#' +#' @seealso \code{\link{arx}}, \code{\link{iv4}} +#' +#' @export +iv <- function(z,order=c(0,1,0),x=NULL){ + if(is.null(x)){ + # Initial Guess using ARX + mod_arx <- arx(z,order) + x <- matrix(sim(mod_arx$sys,inputData(z))) + } + + ivcompute(z,x,order) +} + +#' @import signal +ivcompute <- function(z,x,order,filt=NULL){ + y <- outputData(z); u <- inputData(z); N <- dim(y)[1] + na <- order[1];nb <- order[2]; nk <- order[3] + + nb1 <- nb+nk-1 ; n <- max(na,nb1); df <- N-na-nb + yout <- apply(y,2,padZeros,n=n); + xout <- apply(x,2,padZeros,n=n); + uout <- apply(u,2,padZeros,n=n); + + # Regressors + reg <- function(i) { + if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1 + c(-yout[i-1:na,,drop=T],uout[v,,drop=T]) + } + phi <- t(sapply(n+1:(N+n),reg)) + Y <- yout[n+1:(N+n),,drop=F] + + # Generating IVs + ivx <- function(i) { + if(nk==0) v <- i-0:(nb-1) else v <- i-nk:nb1 + c(-xout[i-1:na,,drop=T],uout[v,,drop=T]) + } + psi <- t(sapply(n+1:(N+n),ivx)) + + l <- list(phi,psi,Y) + if(!is.null(filt)){ + appfilt <- function(x) apply(x,2,signal::filter,filt=filt) + l <- lapply(l, appfilt) + } + phif <- l[[1]]; psif <- l[[2]]; Yf <- l[[3]] + + # Estimator + lhs <- t(psif)%*%phif; lhsinv <- solve(lhs) + theta <- lhsinv%*%t(psif)%*%Yf + + # Residuals + ypred <- (phi%*%theta)[1:N,,drop=F] + e <- y-ypred + sigma2 <- norm(e,"2")^2/df + vcov <- sigma2*solve(t(phi)%*%phi) + + model <- idpoly(A = c(1,theta[1:na]),B = theta[na+1:nb],ioDelay = nk, + Ts=deltat(z),noiseVar = sqrt(sigma2),unit=unit) + + estpoly(sys = model, + stats=list(vcov = vcov, sigma = sqrt(sigma2),df = df), + fitted.values=ypred,residuals=e,call=match.call(),input=u) +} + +#' ARX model estimation using four-stage instrumental variable method +#' +#' Estimates an ARX model of the specified order from input-output data using +#' the instrument variable method. The estimation algorithm is insensitive to +#' the color of the noise term. +#' +#' @param z an idframe object containing the data +#' @param order Specification of the orders: the three integer components +#' (na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) +#' and the input-output delay +#' +#' @details +#' Estimation is performed in 4 stages. The first stage uses the arx function. The resulting model generates the +#' instruments for a second-stage IV estimate. The residuals obtained from this model are modeled using a sufficently +#' high-order AR model. At the fourth stage, the input-output data is filtered through this AR model and then subjected +#' to the IV function with the same instrument filters as in the second stage. +#' +#' @return +#' An object of class \code{estpoly} containing the following elements: +#' \item{sys}{an \code{idpoly} object containing the +#' fitted ARX coefficients} +#' \item{fitted.values}{the predicted response} +#' \item{residuals}{the residuals} +#' \item{input}{the input data used} +#' \item{call}{the matched call} +#' \item{stats}{A list containing the following fields: \cr +#' \code{vcov} - the covariance matrix of the fitted coefficients \cr +#' \code{sigma} - the standard deviation of the innovations\cr +#' \code{df} - the residual degrees of freedom} +#' +#' @references +#' Lennart Ljung (1999), \emph{System Identification: Theory for the User}, +#' 2nd Edition, Prentice Hall, New York. Section 15.3 +#' +#' @examples +#' mod_dgp <- idpoly(A=c(1,-0.5),B=c(0.6,-.2),C=c(1,0.6),ioDelay = 2,noiseVar = 0.1) +#' u <- idinput(400,"prbs") +#' y <- sim(mod_dgp,u,addNoise=TRUE) +#' z <- idframe(y,u) +#' mod_iv4 <- iv4(z,c(1,2,2)) +#' +#' @seealso \code{\link{arx}}, \code{\link{iv4}} +#' @export +iv4 <- function(z,order=c(0,1,0)){ + na <- order[1]; nb <- order[2] + # Steps 1-2 + mod_iv <- iv(z,order) + + # Step 3 + w <- resid(mod_iv) + mod_ar <- ar(w,aic = F,order.max =na+nb) + Lhat <- signal::Arma(b=c(1,-mod_ar$ar),a=1) + + # Step 4 + x2 <- matrix(sim(mod_iv$sys,inputData(z))) + ivcompute(z,x2,order,Lhat) +}
\ No newline at end of file diff --git a/R/nonparam.R b/R/nonparam.R new file mode 100644 index 0000000..8e8e412 --- /dev/null +++ b/R/nonparam.R @@ -0,0 +1,271 @@ +#' Estimate Impulse Response Coefficients +#' +#' \code{impulseest} is used to estimate impulse response coefficients from +#' the data +#' +#' @param x an object of class \code{idframe} +#' @param M Order of the FIR Model (Default:\code{30}) +#' @param K Transport delay in the estimated impulse response +#' (Default:NULL) +#' @param regul Parameter indicating whether regularization should be +#' used. (Default:\code{FALSE}) +#' @param lambda The value of the regularization parameter. Valid only if +#' \code{regul=TRUE}. (Default:\code{1}) +#' +#' @details +#' The IR Coefficients are estimated using linear least squares. Future +#' Versions will provide support for multivariate data. +#' +#' @references +#' Arun K. Tangirala (2015), \emph{Principles of System Identification: +#' Theory and Practice}, CRC Press, Boca Raton. Sections 17.4.11 and 20.2 +#' +#' @seealso \code{\link{step}} +#' +#' @examples +#' uk <- rnorm(1000,1) +#' yk <- filter (uk,c(0.9,-0.4),method="recursive") + rnorm(1000,1) +#' data <- idframe(output=data.frame(yk),input=data.frame(uk)) +#' fit <- impulseest(data) +#' impulseplot(fit) +#' +#' @export +impulseest <- function(x,M=30,K=NULL,regul=F,lambda=1){ + + N <- dim(x$output)[1] + if(is.null(K)) + K <- rep(0,nInputSeries(x)*nOutputSeries(x)) + + out <- rep(list(0),length(K)) + + for(i in seq(nOutputSeries(x))){ + for(j in seq(nInputSeries(x))){ + index <- (i-1)*nInputSeries(x)+j + out[[index]] <- impulsechannel(outputData(x)[,i,drop=F], + inputData(x)[,j,drop=F],N,M, + K[index],regul,lambda) + } + } + out$ninputs <- nInputSeries(x) + out$noutputs <- nOutputSeries(x) + class(out) <- "impulseest" + return(out) +} + +impulsechannel <- function(y,u,N,M,K=0,regul=F,lambda=1){ + ind <- (M+K+1):N + z_reg <- function(i) u[(i-K):(i-M-K),] + Z <- t(sapply(ind,z_reg)) + Y <- y[ind,] + + # Dealing with Regularization + if(regul==F){ + # Fit Linear Model and find standard errors + fit <- lm(Y~Z-1) + coefficients <- coef(fit); residuals <- resid(fit) + } else{ + inner <- t(Z)%*%Z + lambda*diag(dim(Z)[2]) + pinv <- solve(inner)%*% t(Z) + coefficients <- pinv%*%Y + residuals <- Y - Z%*%coefficients + } + df <- nrow(Z)-ncol(Z);sigma2 <- sum(residuals^2)/df + vcov <- sigma2 * solve(t(Z)%*%Z) + se <- sqrt(diag(vcov)) + + out <- list(coefficients=coefficients,residuals=residuals,lags=K:(M+K), + x=colnames(u),y=colnames(y),se = se) + out +} + +#' Impulse Response Plots +#' +#' Plots the estimated IR coefficients along with the significance limits +#' at each lag. +#' +#' @param model an object of class \code{impulseest} +#' @param sd standard deviation of the confidence region (Default: \code{2}) +#' +#' @seealso \code{\link{impulseest}},\code{\link{step}} +#' @import ggplot2 +#' +#' @export +impulseplot <- function(model,sd=2){ + plotseq <- seq(model$noutputs*model$ninputs) + g <- vector("list",model$nin*model$nout) + + for(i in plotseq){ + z <- model[[i]] + lim <- z$se*sd + yindex <- (i-1)%/%model$nin + 1;uindex <- i-model$nin*(yindex-1) + df <- data.frame(x=z$lags,y=coef(z),lim=lim) + g[[i]] <- with(df,ggplot(df,aes(x,y))+ + geom_segment(aes(xend=x,yend=0))+geom_hline(yintercept = 0) + + geom_point(size=2) + ggtitle(paste("From",z$x,"to",z$y))+ + geom_line(aes(y=lim),linetype="dashed",colour="steelblue") + + geom_line(aes(y=-lim),linetype="dashed",colour="steelblue") + + ggplot2::theme_bw(14) + ylab(ifelse(uindex==1,"IR Coefficients","")) + + xlab(ifelse(yindex==model$nout,"Lags","")) + + theme(axis.title=element_text(size=12,color = "black",face = "plain"), + title=element_text(size=9,color = "black",face="bold")) + + scale_x_continuous(expand = c(0.01,0.01))) + } + multiplot(plotlist=g,layout=plotseq) +} + + +#' Step Response Plots +#' +#' Plots the step response of a system, given the IR model +#' +#' @param model an object of class \code{impulseest} +#' +#' @seealso \code{\link{impulseest}} +#' +#' @examples +#' uk <- rnorm(1000,1) +#' yk <- filter (uk,c(0.9,-0.4),method="recursive") + rnorm(1000,1) +#' data <- idframe(output=data.frame(yk),input=data.frame(uk)) +#' fit <- impulseest(data) +#' step(fit) +#' +#' @import ggplot2 +#' @export +step <- function(model){ + plotseq <- seq(model$noutputs*model$ninputs) + g <- vector("list",model$nin*model$nout) + + for(i in plotseq){ + z <- model[[i]] + stepResp <- cumsum(coef(z)) + yindex <- (i-1)%/%model$nin + 1;uindex <- i-model$nin*(yindex-1) + df <- data.frame(x=z$lags,y=stepResp) + g[[i]] <- with(df,ggplot(df,aes(x,y))+ + geom_step() + ggtitle(paste("From",z$x,"to",z$y)) + + theme_bw(14) + ylab(ifelse(uindex==1,"Step Response","")) + + xlab(ifelse(yindex==model$nout,"Lags","")) + + theme(axis.title=element_text(size=12,color = "black",face = "plain"), + title=element_text(size=9,,color = "black",face="bold"))) + } + multiplot(plotlist=g,layout=plotseq) +} + +#' Estimate frequency response +#' +#' Estimates frequency response and noise spectrum from data with +#' fixed resolution using spectral analysis +#' +#' @param x an \code{idframe} object +#' @param winsize lag size of the Hanning window (Default: \code{min +#' (length(x)/10,30)}) +#' @param freq frequency points at which the response is evaluated +#' (Default: \code{seq(1,128)/128*pi/Ts}) +#' +#' @return +#' an \code{idfrd} object containing the estimated frequency response +#' and the noise spectrum +#' +#' @references +#' Arun K. Tangirala (2015), \emph{Principles of System Identification: +#' Theory and Practice}, CRC Press, Boca Raton. Sections 16.5 and 20.4 +#' +#' @examples +#' data(arxsim) +#' frf <- spa(arxsim) +#' +#' @import signal +#' @export +spa <- function(x,winsize=NULL,freq=NULL){ + N <- dim(x$out)[1] + nout <- nOutputSeries(x); nin <- nInputSeries(x) + + if(is.null(winsize)) winsize <- min(N/10,30) + if(is.null(freq)) freq <- (1:128)/128*pi/deltat(x) + M <- winsize + + Ryu <- mult_ccf(x$out,x$input,lag.max = M) + Ruu <- mult_ccf(x$input,x$input,lag.max=M) + Ryy <- mult_ccf(x$out,x$out,lag.max = M) + + cov2spec <- function(omega,R,M){ + seq1 <- exp(-1i*(-M:M)*omega) + sum(R*signal::hanning(2*M+1)*seq1) + } + + G <- array(0,c(nout,nin,length(freq))) + spec <- array(0,c(nout,1,length(freq))) + for(i in 1:nout){ + phi_y <- sapply(freq,cov2spec,Ryy[i,i,],M) + temp <- phi_y + for(j in 1:nin){ + phi_yu <- sapply(freq,cov2spec,Ryu[i,j,],M) + phi_u <- sapply(freq,cov2spec,Ruu[j,j,],M) + G[i,j,] <- phi_yu/phi_u + temp <- temp - phi_yu*Conj(phi_yu)/phi_u + } + spec[i,1,] <- temp + } + out <- idfrd(G,matrix(freq),deltat(x),spec) + return(out) +} + +mult_ccf <- function(X,Y=NULL,lag.max=30){ + N <- dim(X)[1]; nx <- dim(X)[2] + ny <- ifelse(is.null(Y),nx,dim(Y)[2]) + + ccvfij <- function(i,j,lag=30) ccf(X[,i],Y[,j],plot=F,lag.max =lag, + type="covariance") + Xindex <- matrix(sapply(1:nx,rep,nx),ncol=1)[,1] + temp <- mapply(ccvfij,i=Xindex,j=rep(1:ny,ny), + MoreArgs = list(lag=lag.max)) + + ccfextract <- function(i,l) l[,i]$acf + temp2 <- t(sapply(1:(nx*ny),ccfextract,l=temp)) + dim(temp2) <- c(nx,ny,2*lag.max+1) + return(temp2) +} + +#' Estimate empirical transfer function +#' +#' Estimates the emperical transfer function from the data by taking the +#' ratio of the fourier transforms of the output and the input variables +#' +#' @param data an object of class \code{idframe} +#' @param n frequency spacing (Default: \code{128}) +#' +#' @return +#' an \code{idfrd} object containing the estimated frequency response +#' +#' @references +#' Arun K. Tangirala (2015), \emph{Principles of System Identification: +#' Theory and Practice}, CRC Press, Boca Raton. Sections 5.3 and 20.4.2 +#' +#' @seealso \code{\link[stats]{fft}} +#' +#' @examples +#' data(arxsim) +#' frf <- etfe(arxsim) +#' +#' @export +etfe <- function(data,n=128){ + y <- data$output + u <- data$input + N <- dim(data$output)[1] + if(N < n){ + n=N + } + v=seq(1,N,length.out = n) + temp <- cbind(data$output[v,],data$input[v,]) + tempfft <- mvfft(temp)/dim(temp)[1] + G <- comdiv(tempfft[,1],tempfft[,2]) + resp = G[1:ceiling(length(G)/2)] + frequency <- matrix(seq( 1 , ceiling(n/2) ) * pi / floor(n/2) / deltat(data)) + out <- idfrd(respData = resp,freq=frequency,Ts=data$Ts) + return(out) +} +comdiv <- function(z1,z2){ + mag1 <- Mod(z1);mag2 <- Mod(z2) + phi1 <- Arg(z1); phi2 <- Arg(z2) + + complex(modulus=mag1/mag2,argument=signal::unwrap(phi1-phi2)) +}
\ No newline at end of file diff --git a/R/poly.R b/R/poly.R new file mode 100644 index 0000000..528a599 --- /dev/null +++ b/R/poly.R @@ -0,0 +1,196 @@ +#' Polynomial model with identifiable parameters +#' +#' Creates a polynomial model with identifiable coefficients +#' +#' @param A autoregressive coefficients +#' @param B,F1 coefficients of the numerator and denominator respectively +#' of the deterministic model between the input and output +#' @param C,D coefficients of the numerator and denominator respectively +#' of the stochastic model +#' @param ioDelay the delay in the input-output channel +#' @param Ts sampling interval +#' @param noiseVar variance of the white noise source (Default=\code{1}) +#' @param intNoise Logical variable indicating presence or absence of integrator +#' in the noise channel (Default=\code{FALSE}) +#' @param unit time unit (Default=\code{"seconds"}) +#' +#' @details +#' Discrete-time polynomials are of the form +#' \deqn{ +#' A(q^{-1}) y[k] = \frac{B(q^{-1})}{F1(q^{-1})} u[k] + +#' \frac{C(q^{-1})}{D(q^{-1})} e[k] +#' } +#' +#' @examples +#' # define output-error model +#' mod_oe <- idpoly(B=c(0.6,-0.2),F1=c(1,-0.5),ioDelay = 2,Ts=0.1, +#' noiseVar = 0.1) +#' +#' # define box-jenkins model with unit variance +#' B <- c(0.6,-0.2) +#' C <- c(1,-0.3) +#' D <- c(1,1.5,0.7) +#' F1 <- c(1,-0.5) +#' mod_bj <- idpoly(1,B,C,D,F1,ioDelay=1) +#' +#' @export +idpoly <- function(A=1,B=1,C=1,D=1,F1=1,ioDelay=0,Ts=1, + noiseVar=1,intNoise = F,unit = c("seconds","minutes", + "hours","days")[1]){ + Bindex <- which.max(B!=0) + ioDelay <- ifelse(Bindex-1==0,ioDelay,Bindex-1) + B <- B[Bindex:length(B)] + out <- list(A= A,B=B,C=C,D=D,F1=F1,ioDelay = ioDelay,Ts=Ts, + noiseVar=noiseVar,unit=unit,intNoise = intNoise) + out$type <- typecheck(out) + class(out) <- "idpoly" + return(out) +} + +typecheck <- function(x){ + y <- lapply(x[1:5],checkUnity) + if(y$A){ + out <- if(y$C && y$D) if(y$F1) "fir" else "oe" else "bj" + } else{ + if(y$D && y$F1){ + if(x$intNoise){ + out <- if(y$C) "ari" else "arima" + } else{ + out <- if(y$C) "ar" else "arma" + } + if(!y$B) out <- paste(out,"x",sep="") + } else{ + out <- "polynomial" + } + } + out +} + +checkUnity <- function(x){ + out <- if(length(x)==1 && x==1) TRUE else FALSE +} + +#' @export +print.idpoly <- function(x,se=NULL,dig=3,...){ + main <- paste("Discrete-time",toupper(x$type),"model:") + if(x$type=="oe" || x$type=="bj"){ + main <- paste(main,"y[k] = B(z)/F(z) u[k] +") + if(x$type=="bj"){ + polyExp <- if(!checkUnity(x$C)) "C(z)/D(z)" else "1/D(z)" + if(x$intNoise==T) polyExp <- paste(polyExp,"1/(1-z^{-1})") + main <- paste(main,polyExp,"e[k] \n\n") + } + } else{ + main <- paste(main,"A(z)y[k] =") + if(!checkUnity(x$B)) main <- paste(main,"B(z)u[k] +") + if(checkUnity(x$C)){ + Cexp <- if(!x$intNoise) "" else "1/(1-z^{-1})" + }else{ + Cexp <- "C(z)" + if(x$intNoise) Cexp <- paste(Cexp,"/(1-z^{-1})",sep="") + } + main <- paste(main,Cexp,"e[k] \n\n") + } + cat(main) + + # Printing Standard error sequence + j=1 + print_se <- function(se){ + if(!is.null(se)){ + cat(" (+/- ",round(se[j],dig),") ",sep = "") + j <<- j+1 + } + } + + if(length(x$A)>1){ + cat("A(q^{-1}) = ") + for(i in seq_along(x$A)){ + if(i-1==0){ + cat(round(x$A[i],dig)) + } else{ + if(x$A[i]>0) cat(" + ") else cat("- ") + + if(!(abs(x$A[i])==1)) cat(abs(round(x$A[i],dig))) + print_se(se) + cat("q^{-",i-1,"}",sep="") + } + cat("\t") + } + cat("\n") + } + + cat("B(q^{-1}) = ") + for(i in seq_along(x$B)){ + if(i+x$ioDelay-1==0){ + cat(round(x$B[i],dig)) + } else{ + + if(!((x$ioDelay!=0) && (i==1))){ + if(x$B[i]>0) cat(" + ") else cat("- ") + } else{ + if(x$B[i]<0) cat("-") + } + + if(!(abs(x$B[i])==1)) cat(abs(round(x$B[i],dig))) + print_se(se) + cat("q^{-",i+x$ioDelay-1,"}",sep="") + } + cat("\t") + } + cat("\n") + + if(length(x$C)>1){ + cat("C(q^{-1}) = ") + for(i in seq_along(x$C)){ + if(i-1==0){ + cat(round(x$C[i],dig)) + } else{ + if(x$C[i]>0) cat(" + ") else cat("- ") + + if(!(abs(x$C[i])==1)) cat(abs(round(x$C[i],dig))) + print_se(se) + cat("q^{-",i-1,"}",sep="") + } + cat("\t") + } + cat("\n") + } + + if(length(x$D)>1){ + cat("D(q^{-1}) = ") + for(i in seq_along(x$D)){ + if(i-1==0){ + cat(round(x$D[i],dig)) + } else{ + if(x$D[i]>0) cat(" + ") else cat("- ") + + if(!(abs(x$D[i])==1)) cat(abs(round(x$D[i],dig))) + print_se(se) + cat("q^{-",i-1,"}",sep="") + } + cat("\t") + } + cat("\n") + } + + if(length(x$F1)>1){ + cat("F(q^{-1}) = ") + for(i in seq_along(x$F1)){ + if(i-1==0){ + cat(round(x$F1[i],dig)) + } else{ + if(x$F1[i]>0) cat(" + ") else cat("- ") + + if(!(abs(x$F1[i])==1)) cat(abs(round(x$F1[i],dig))) + print_se(se) + cat("q^{-",i-1,"}",sep="") + } + cat("\t") + } + } + cat("\n") +} + +params <- function(x){ + c(x$A[-1],x$B,x$C[-1],x$D[-1],x$F1[-1]) +}
\ No newline at end of file diff --git a/R/predict.R b/R/predict.R new file mode 100644 index 0000000..b56fa8f --- /dev/null +++ b/R/predict.R @@ -0,0 +1,120 @@ +predict.idpoly <- function(object,data,nahead=1){ + y <- outputData(data); u<- inputData(data) + G <- signal::Arma(b=c(rep(0,object$ioDelay),object$B), + a= as.numeric(polynom::polynomial(object$A)* + polynom::polynomial(object$F1))) + det_sys <- as.numeric(signal::filter(G,u)) + if(object$type=="oe" || nahead==Inf){ + ypred <- det_sys + } else{ + Hden <- as.numeric(polynom::polynomial(object$A)*polynom::polynomial(object$D)) + Hinv <- signal::Arma(b=Hden,a=object$C) + filtered <- as.numeric(signal::filter(Hinv,as.numeric(y)-det_sys)) + if(nahead!=1){ + H <- as.numeric(polynom::polynomial(object$C)*polyinv(Hden,nahead)) + Hl <- signal::Ma(H[1:nahead]) + filtered <- as.numeric(signal::filter(Hl,filtered)) + } + ypred <- as.numeric(y) - filtered + } + + matrix(ts(ypred,start=start(data),deltat=deltat(data))) +} + +polyinv <- function(x,k){ + gamma <- 1/Re(polyroot(x)) + + inverse <- function(y,k){ + sapply(1:k-1,function(i) y^i) + } + z <- lapply(lapply(gamma,inverse,k=2),polynom::polynomial) + temp = z[[1]] + if(length(z)>1){ + for(i in 2:length(z)){ + temp = temp*z[[i]] + } + } + temp +} + +#' Predictions of identified model +#' +#' Predicts the output of an identified model (\code{estpoly}) object K steps ahead. +#' +#' @param object \code{estpoly} object containing the identified model +#' @param newdata optional dataset to be used for predictions. If not supplied, +#' predictions are made on the training set. +#' @param nahead number of steps ahead at which to predict (Default:1). For infinite- +#' step ahead predictions or pure simulation, supply \code{Inf}. +#' @param \ldots other arguments +#' +#' @return +#' Time-series containing the predictions +#' +#' @examples +#' data(arxsim) +#' mod1 <- oe(arxsim,c(2,1,1)) +#' Yhat <- predict(mod1,arxsim) # 1-step ahead predictions +#' Yhat_2 <- predict(mod1,arxsim,nahead=2) # 2-step ahead predictions +#' Yhat_inf <- predict(mod1,arxsim,nahead=Inf) # Infinite-step ahead predictions +#' +#' @references +#' Arun K. Tangirala (2015), \emph{Principles of System Identification: Theory +#' and Practice}, CRC Press, Boca Raton. Chapter 18 +#' +#' @export +predict.estpoly <- function(object,newdata=NULL,nahead=1,...){ + if(is.null(newdata)&& nahead==1){ + return(matrix(fitted(object))) + } else{ + model <- object$sys + if(is.null(newdata)){ + y <- fitted(object)+resid(object) + u <- object$input + z <- idframe(y,u,Ts = deltat(y),start=start(y)) + } else{ + z <- newdata + } + predict(model,z,nahead) + } +} + +#' Compare the measured output and the predicted output(s) +#' +#' Plots the output predictions of model(s) superimposed over validation data, +#' data, for comparison. +#' +#' @param data validation data in the form of an \code{idframe} object +#' @param nahead number of steps ahead at which to predict (Default:1). For infinite- +#' step ahead predictions, supply \code{Inf}. +#' @param \ldots models whose predictions are to be compared +#' +#' @examples +#' data(arxsim) +#' mod1 <- arx(arxsim,c(1,2,2)) +#' mod2 <- oe(arxsim,c(2,1,1)) +#' compare(arxsim,nahead=Inf,mod1,mod2) +#' +#' @seealso \code{\link{predict.estpoly}} for obtaining model predictions +#' @import ggplot2 reshape2 +#' @export +compare <- function(data,nahead=1,...){ + loadNamespace("ggplot2") + # Input Validation + input_list <- as.list(substitute(list(...)))[-1] + dots <- list(...) + if(is.null(dots)) stop("No model supplied") + + Y <- sapply(dots,predict,newdata=data,nahead=nahead) + nrmse <- sapply(dots,FUN = function(x) fitch(x)$FitPer) + temp <- paste(as.character(input_list),paste(round(nrmse,2),"%",sep=""), + sep=": ") + df <- data.frame(Time = as.numeric(time(data)), + Actual=as.numeric(outputData(data)[,1]),Y) + meltdf <- reshape2::melt(df,id="Time") + + with(meltdf,ggplot(meltdf,aes(x=Time,y=value,color=variable,group=variable))+geom_line(size=1)+ + ggtitle(paste(nahead,"step(s) ahead model predictions"))+ + theme_bw()+ylab(outputNames(data)) + labs(colour="") + + scale_colour_hue(labels=c("Measured",temp),l=50)) +}
\ No newline at end of file diff --git a/R/preprocess.R b/R/preprocess.R new file mode 100644 index 0000000..b5ba991 --- /dev/null +++ b/R/preprocess.R @@ -0,0 +1,192 @@ +#' Remove offsets and linear trends +#' +#' Removes offsets or trends from data +#' @param x an object of class \code{idframe} +#' @param type argument indicating the type of trend to be removed (Default=\code{0}) +#' \itemize{ +#' \item type=\code{0}: Subtracts mean value from each signal +#' \item type=\code{1}: Subtracts a linear trend (least-squres fit) +#' \item type=\code{trInfo} object: Subtracts a trend specified by the object +#' } +#' +#' @return +#' A list containing two objects: the detrended data and the trend information +#' +#' @details +#' +#' \code{R} by default doesn't allow return of multiple objects. The \code{\%=\%} +#' operator and \code{g} function in this package facillitate this behaviour. See +#' the examples section for more information. +#' +#' @aliases trInfo +#' +#' @examples +#' data(cstr) +#' datatrain <- dataSlice(cstr,end=4500) +#' datatest <- dataSlice(cstr,4501) +#' g(Ztrain,tr) %=% detrend(datatrain) # Remove means +#' g(Ztest) %=% detrend(datatest,tr) +#' +#' @seealso \code{\link[stats]{lm}} +#' @export +detrend <- function(x,type=0){ + z <- x + reg <- time(x) + if(class(type)=="trInfo"){ # remove custom trend + if(nOutputSeries(x)!=0){ + fit <- sweep(sweep(matrix(rep(reg,nOutputSeries(x)),ncol=nOutputSeries(x)), + 2,type$OutputSlope,"*"),2,type$OutputOffset,"+") + z$output <- x$output-fit + } + if(nInputSeries(x)!=0){ + fit <- sweep(sweep(matrix(rep(reg,nInputSeries(x)),ncol=nInputSeries(x)), + 2,type$InputSlope,"*"),2,type$InputOffset,"+") + z$input <- x$input-fit + } + tinfo <- type + } else if(type == 0){ # remove means + tinfo <- trInfo() + if(nOutputSeries(x)!=0){ + outputData(z) <- apply(outputData(x),2,scale,T,F) + tinfo$OutputOffset <- colMeans(x$output) + tinfo$OutputSlope <- rep(0,nOutputSeries(x)) + } + + if(nInputSeries(x)!=0){ + inputData(z) <- apply(inputData(x),2,scale,T,F) + tinfo$InputOffset <- colMeans(x$input) + tinfo$InputSlope <- rep(0,nInputSeries(x)) + } + } else if(type==1){ + formula <- X ~ reg + + # Function which performs linear regression across every column + multilm <- function(x,formula,time){ + l <- lapply(as.list(x),function(x) data.frame(X=x,reg=time)) + trend <- lapply(l,function(x) lm(formula,data=x)) + trend + } + + tinfo <- trInfo() + if(nOutputSeries(x)!=0){ + output_trend <- multilm(outputData(x),formula,reg) + outputData(z) <- ts(sapply(output_trend,resid),start=reg[1], + end=tail(reg,n=1),deltat=deltat(x)) + out_coefs <- sapply(output_trend,coef) + tinfo$OutputOffset <- out_coefs[1,] + tinfo$OutputSlope <- out_coefs[2,] + } + + if(nInputSeries(x)!=0){ + input_trend <- multilm(inputData(x),formula,reg) + inputData(z) <- ts(sapply(input_trend,resid),start=reg[1], + end=tail(reg,n=1),deltat=deltat(x)) + in_coefs <- sapply(input_trend,coef) + tinfo$InputOffset <- in_coefs[1,] + tinfo$InputSlope <- in_coefs[2,] + } + } else{ + stop("Error: Invalid trend type") + } + list(z,tinfo) +} + +#' @export +trInfo <- function(InputOffset=numeric(0),OutputOffset=numeric(0), + InputSlope=numeric(0),OutputSlope=numeric(0)){ + l <- list(InputOffset=InputOffset,OutputOffset=OutputOffset, + InputSlope=InputSlope,OutputSlope=OutputSlope) + class(l) <- "trInfo" + l +} + +#' Replace Missing Data by Interpolation +#' +#' Function for replacing missing values with interpolated ones. This is an +#' extension of the \code{na.approx} function from the \code{zoo} package. +#' The missing data is indicated using the value \emph{NA}. +#' +#' @param data an object of class \code{idframe} +#' @return +#' data (an idframe object) with missing data replaced. +#' +#' @seealso \code{\link[zoo]{na.approx}} +#' +#' @examples +#' data(cstr_mis) +#' summary(cstr_mis) # finding out the number of NAs +#' cstr <- misdata(cstr_mis) +#' +#' @importFrom zoo na.approx +#' @export +misdata <- function(data){ + if (!requireNamespace("zoo", quietly = TRUE)) { + stop("Package zoo needed for this function to work. Please install it.", + call. = FALSE) + } + + f <- function(var,start,end,Ts){ + time_range <- range(time(var)) + start <- time_range[1];end <- time_range[2] + Ts <- stats::deltat(var) + var <- ts(data=var,start=start,end=end,deltat=Ts) + out <- zoo::na.approx(var,na.rm=F) + return(as.numeric(out)) + } + + Z <- data + outputData(Z) <- apply(outputData(data),2,f) + inputData(Z) <- apply(inputData(data),2,f) + Z +} + + +#' Subset or Resample idframe data +#' +#' \code{dataSlice} is a subsetting method for objects of class \code{idframe}. It +#' extracts the subset of the object \code{data} observed between indices \code{start} +#' and \code{end}. If a frequency is specified, the series is then re-sampled at the +#' new frequency. +#' +#' @param data an object of class \code{idframe} +#' @param start the start index +#' @param end the end index +#' @param freq fraction of the original frequency at which the series +#' to be sampled. +#' +#' @details +#' The dataSlice function extends the \code{\link[stats]{window}} +#' function for idframe objects +#' +#' @return an idframe object +#' +#' @examples +#' data(cstr) +#' cstrsub <- dataSlice(cstr,start=200,end=400) # extract between indices 200 and 400 +#' cstrTrain <- dataSlice(cstr,end=4500) # extract upto index 4500 +#' cstrTest <- dataSlice(cstr,start=6501) # extract from index 6501 till the end +#' cstr_new <- dataSlice(cstr,freq=0.5) # resample data at half the original frequency +#' +#' @seealso \code{\link[stats]{window}} +#' @export +dataSlice <- function(data,start=NULL,end=NULL,freq=NULL){ + # check if the class is correct + if(class(data)!='idframe') + stop("Not an idframe data") + + indexWindow <- function(y,start,end,freq){ + Y <- matrix(y,ncol=ncol(y)); z <- as.vector(time(y)) + Y <- window(Y,start=start,end=end,frequency=freq) + zw <- window(z,start=start,end=end,frequency=freq) + temp <- ts(Y,start=zw[1],end=tail(zw,n=1),deltat=diff(zw)[1]) + colnames(temp) <- colnames(y) + temp + } + if(nOutputSeries(data)!=0) + outputData(data) <- indexWindow(outputData(data),start,end,freq) + + if(nInputSeries(data)!=0) + inputData(data) <- indexWindow(inputData(data),start,end,freq) + + return(data) +}
\ No newline at end of file diff --git a/R/rarx.R b/R/rarx.R new file mode 100644 index 0000000..7f48888 --- /dev/null +++ b/R/rarx.R @@ -0,0 +1,79 @@ +#' Estimate parameters of ARX recursively +#' +#' Estimates the parameters of a single-output ARX model of the +#' specified order from data using the recursive weighted least-squares +#' algorithm. +#' +#' @param x an object of class \code{idframe} +#' @param order Specification of the orders: the three integer components +#' (na,nb,nk) are the order of polynolnomial A, (order of polynomial B + 1) and +#' the input-output delay +#' @param lambda Forgetting factor(Default=\code{0.95}) +#' +#' @return +#' A list containing the following objects +#' \describe{ +#' \item{theta}{Estimated parameters of the model. The \eqn{k^{th}} +#' row contains the parameters associated with the \eqn{k^{th}} +#' sample. Each row in \code{theta} has the following format: \cr +#' theta[i,:]=[a1,a2,...,ana,b1,...bnb] +#' } +#' \item{yhat}{Predicted value of the output, according to the +#' current model - parameters based on all past data} +#' } +#' +#' @references +#' Arun K. Tangirala (2015), \emph{Principles of System Identification: +#' Theory and Practice}, CRC Press, Boca Raton. Section 25.1.3 +#' +#' Lennart Ljung (1999), \emph{System Identification: Theory for the User}, +#' 2nd Edition, Prentice Hall, New York. Section 11.2 +#' @examples +#' Gp1 <- idpoly(c(1,-0.9,0.2),2,ioDelay=2,noiseVar = 0.1) +#' Gp2 <- idpoly(c(1,-1.2,0.35),2.5,ioDelay=2,noiseVar = 0.1) +#' uk = idinput(2044,'prbs',c(0,1/4)); N = length(uk); +#' N1 = round(0.35*N); N2 = round(0.4*N); N3 = N-N1-N2; +#' yk1 <- sim(Gp1,uk[1:N1],addNoise = TRUE) +#' yk2 <- sim(Gp2,uk[N1+1:N2],addNoise = TRUE) +#' yk3 <- sim(Gp1,uk[N1+N2+1:N3],addNoise = TRUE) +#' yk <- c(yk1,yk2,yk3) +#' z <- idframe(yk,uk,1) +#' g(theta,yhat) %=% rarx(z,c(2,1,2)) +#' +#' @export +rarx <- function(x,order=c(1,1,1),lambda=0.95){ + y <- outputData(x); u <- inputData(x) + N <- dim(y)[1] + na <- order[1];nb <- order[2]; nk <- order[3] + nb1 <- nb+nk-1 ; n <- max(na,nb1); df <- N-na-nb + + yout <- apply(y,2,padZeros,n=n) + uout <- apply(u,2,padZeros,n=n) + + uindex <- nk:nb1 + if(na!=0) yindex <- 1:na + + reg <- function(i) { + # regressor + temp <- numeric(0) + if(na!=0) temp <- c(temp,-yout[i-yindex,]) + phi <- c(temp,uout[i-uindex,]) + phi + } + + # R0 <- reg(n+1)%*%t(reg(n+1)) + # Plast <- solve(R0) + Plast <- 10^4*diag(na+nb) + theta <- matrix(0,N+1,na+nb) + yhat <- y + + for(i in 1:N){ + temp <- reg(n+i) + yhat[i,] <- t(temp)%*%t(theta[i,,drop=FALSE]) + eps_i <- y[i,,drop=FALSE] - yhat[i,,drop=FALSE] + kappa_i <- Plast%*%temp/(lambda+t(temp)%*%Plast%*%temp)[1] + theta[i+1,] <- t(t(theta[i,,drop=F])+eps_i[1]*kappa_i) + Plast <- (diag(na+nb)-kappa_i%*%t(temp))%*%Plast/lambda + } + list(theta=theta[1+1:N,],yhat=yhat) +}
\ No newline at end of file diff --git a/R/readData.R b/R/readData.R new file mode 100644 index 0000000..33c1dcb --- /dev/null +++ b/R/readData.R @@ -0,0 +1,73 @@ +#' Data input into a idframe object +#' +#' Read the contents of a data.frame/matrix into a \code{idframe} object. +#' +#' @param data a \code{data.frame} object +#' @param ninputs the number of input columns. (Default: 0) +#' @param Ts sampling interval (Default: 1) +#' @param unit Time Unit (Default: "seconds") +#' +#' +#' @return an idframe object +#' @examples +#' data(cstrData) +#' data <- read.idframe(cstrData,ninputs=1,Ts= 1,unit="minutes") +#' +#' @export +read.idframe <- function(data,ninputs=NULL,Ts = 1, + unit=c("seconds","minutes","hours", + "days")[1]){ + outIndex <- 1:dim(data)[2]; inputs <- NULL + if(ninputs!=0){ + inputs <- data[,1:ninputs,drop=F] + outIndex <- seq(ninputs+1,dim(data)[2],by=1) + } + outputs <- data[,outIndex,drop=F] + + out <- idframe(output=outputs,input=inputs,Ts=Ts,unit=unit) + + return(out) +} + +#' Read the contents of a table-formatted file +#' +#' Read the contents of an file in table format into a \code{idframe} object. +#' +#' @param file the path to the file to read +#' @param sep the field separator character. Values on each line of the file are +#' separated by this character. (Default: \code{","}) +#' @param header a logical value indicating whether the first row corresponding to +#' the first element of the rowIndex vector contains the names of the variables. +#' (Default: \code{TRUE}) +#' @param ninputs the number of input columns. (Default: 0) +#' @param Ts sampling interval (Default: 1) +#' @param unit Time Unit (Default: "seconds") +#' @param ... additional arguments to be passed to the \code{\link[utils]{read.table}} function +#' +#' @details +#' +#' The \code{read.table.idframe} function uses the \code{\link[utils]{read.table}} function, +#' provided by the \pkg{utils} package, to read data from a table-formatted file and then calls the +#' \code{\link{read.idframe}} function to read the data into a idframe object +#' +#' @return an idframe object +#' @examples +#' dataMatrix <- data.frame(matrix(rnorm(1000),ncol=5)) +#' colnames(dataMatrix) <- c("u1","u2","y1","y2","y3") +#' write.csv(dataMatrix,file="test.csv",row.names=FALSE) +#' +#' data <- read.table.idframe("test.csv",ninputs=2,unit="minutes") +#' +#' @seealso \code{\link[utils]{read.table}} +#' @export +read.table.idframe <- function(file,header=TRUE,sep=",",ninputs=0, + Ts = 1,unit=c("seconds","minutes","hours", + "days")[1],...){ + + # Read from file (default: csv file) + dat <- read.table(file=file,header=header,sep=sep,...) + + # read from dataframe and return idframe object + out <- read.idframe(dat,ninputs=ninputs,Ts = Ts,unit=unit) + return(out) +}
\ No newline at end of file @@ -0,0 +1,87 @@ +#' Simulate response of dynamic system +#' +#' Simulate the response of a system to a given input +#' +#' @param model the linear system to simulate +#' @param input a vector/matrix containing the input +#' @param addNoise logical variable indicating whether to add noise to the +#' response model. (Default: \code{FALSE}) +#' @param innov an optional times series of innovations. If not supplied (specified +#' as \code{NULL}), gaussian white noise is generated, with the variance specified in +#' the model (Property: \code{noiseVar}) +#' @param seed integer indicating the seed value of the random number generator. +#' Useful for reproducibility purposes. +#' +#' @return +#' a vector containing the simulated output +#' +#' @details +#' The routine is currently built only for SISO systems. Future versions will +#' include support for MIMO systems. +#' +#' @examples +#' # ARX Model +#' u <- idinput(300,"rgs") +#' model <- idpoly(A=c(1,-1.5,0.7),B=c(0.8,-0.25),ioDelay=1, +#' noiseVar=0.1) +#' y <- sim(model,u,addNoise=TRUE) +#' +#' @export +sim <- function(model,input,addNoise = F,innov=NULL,seed=NULL) UseMethod("sim") + +#' @export +sim.default <- function(model,input,addNoise = F,innov=NULL,seed=NULL){ + print("The sim method is not developed for the current class of the object") +} + +#' @import signal polynom +#' @export +sim.idpoly <- function(model,input,addNoise = F,innov=NULL,seed=NULL){ + B <- c(rep(0,model$ioDelay),model$B) + Gden <- as.numeric(polynomial(model$A)*polynomial(model$F1)) + G <- signal::Arma(b=B,a=Gden) + ufk <- signal::filter(G,input) + yk <- as.numeric(ufk) + + if(addNoise){ + n <- ifelse(is.numeric(input),length(input),dim(input)[1]) + if(!is.null(innov)){ + ek <- innov + } else{ + if(!is.null(seed)) set.seed(seed) + ek <- rnorm(n,sd=sqrt(model$noiseVar)) + } + den1 <- as.numeric(polynomial(model$A)*polynomial(model$D)) + filt1 <- Arma(b=model$C,a=den1) + vk <- signal::filter(filt1,ek) + + yk <- yk + as.numeric(vk) + } + + return(yk) +} + +# sim_arx <- function(model,input,ek){ +# na <- length(model$A) - 1; nk <- model$ioDelay; +# nb <- length(model$B) - 1; nb1 <- nb+nk +# n <- max(na,nb1) +# coef <- matrix(c(model$A[-1],model$B),nrow=na+(nb+1)) +# +# if(class(input)=="idframe"){ +# uk <- input$input[,1,drop=T] +# } else if(class(input) %in% c("matrix","data.frame")){ +# uk <- input[,1,drop=T] +# } else if(is.numeric(input)){ +# uk <- input +# } +# +# y <- rep(0,length(input)+n) +# u <- c(rep(0,n),uk) +# +# for(i in n+1:length(uk)){ +# if(nk==0) v <- u[i-0:nb] else v <- u[i-nk:nb1] +# reg <- matrix(c(-(y[i-1:na]),v),ncol=na+(nb+1)) +# y[i] <- reg%*%coef + ek[i] +# } +# return(y[n+1:length(uk)]) +# }
\ No newline at end of file diff --git a/R/sysid.R b/R/sysid.R new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/R/sysid.R diff --git a/R/util.R b/R/util.R new file mode 100644 index 0000000..4667f9a --- /dev/null +++ b/R/util.R @@ -0,0 +1,53 @@ +#' Multiple assignment operator +#' +#' Assign multiple variables from a list or function return object +#' +#' @param l the variables to be assigned +#' @param r the list or function-return object +#' +#' @aliases g +#' +#' @export +'%=%' = function(l,r) UseMethod('%=%') + +# Binary Operator +#' @export +'%=%.lbunch' = function(l,r) { + Envir = as.environment(-1) + +# if (length(r) > length(l)) +# warning("RHS has more args than LHS. Only first", length(l), "used.") + + if (length(l) > length(r)) { + warning("LHS has more args than RHS. RHS will be repeated.") + r <- extendToMatch(r, l) + } + + for (II in 1:length(l)) { + do.call('<-', list(l[[II]], r[[II]]), envir=Envir) + } +} + +# Used if LHS is larger than RHS +extendToMatch <- function(source, destin) { + s <- length(source) + d <- length(destin) + + # Assume that destin is a length when it is a single number and source is not + if(d==1 && s>1 && !is.null(as.numeric(destin))) + d <- destin + + dif <- d - s + if (dif > 0) { + source <- rep(source, ceiling(d/s))[1:d] + } + return (source) +} + +# Grouping the left hand side +#' @export +g = function(...) { + List = as.list(substitute(list(...)))[-1L] + class(List) = 'lbunch' + return(List) +}
\ No newline at end of file |