1 Overview

A dynamic system is said to be “controllable” if the internal, underlying values of the system can be driven or directed toward a specific target level within a finite time interval regardless of its initial values through the use of one or more exogenous inputs. Examples may include the use of insulin (the “input”) to drive the glucose level of a diabetic patient toward a target glucose level. This tutorial implements a type of control theory algorithm, called Linear-Quadratic Controller, using intermediate output from the dynr package. This tutorial is also featured as an example in Gates, Chow, and Molenaar (in progress).

Linear-Quadratic Controller is one special case of a class of control strategies called receding horizon controls (RHCs) or model predictive controls. This special case assumes that the true values of a system’s latent underlying (i.e., state) values are known, but instead, have to be estimated from noisy observed variables. Details can be found in Kwon and Ham (2005), and Molenaar (2010).

References:

Gates, K., Molenaar, P. C. M., & Chow, S.-M. (in progress). Analysis of Intra-individual Variation. New York: Taylor & Francis.

Goodwin, G., Seron, M. M., & de Don, J. A. (2005). Constrained control and estimation: An optimisation approach (1st ed.). London, United Kingdom: Springer-Verlag London Ltd.

Kwon, W. H., & Han, S. H. (2005). Receding horizon control: model predictive control for state models. London: Springer.

Molenaar, P. C. (2010). Note on optimization of individual psychotherapeutic processes. Journal of Mathematical Psychology, 54 (1), 208{213.

Wang, Q., Molenaar, P., Harsh, S., Freeman, K., Xie, J., Gold, C., . . . Ulbrecht, J. (2014a). Personalized state-space modeling of glucose dynamics for type-1 diabetes using continuously monitored glucose, insulin dose, and meal intake: An extended kalman lter approach. Journal of Diabetes Science and Technology, 8 (2), 331-345. Retrieved from http://dst.sagepub.com/ content/8/2/331.abstract doi: 10.1177/1932296814524080

2 Load dependent R packages and set up options

## Loading required package: mvtnorm
## Loading required package: ggplot2
## Loading required package: dynr
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## Loading required package: gridExtra
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
## Loading required package: zeallot

3 Specify function to generate data using a user-specified linear state-space model.

GenerateData = function(T,npad,np,ne,ny,r,psi,lambda,theta,beta,
                        eta_intercept,G,a0,P0,zt,
                        fileObs=NULL, fileState=NULL,shocks=NULL){
yall <- matrix(0,nrow = T*np, ncol = ny)
etaall <- matrix(0,nrow = T*np, ncol = ne)
if (is.null(shocks)){
  shocks = matrix(0,np*T,ne)
}
for (p in 1:np){
  etaC      <- matrix(0, nrow = ne, ncol = T + npad) # Set up matrix for contemporaneous variables. 
  etaC[,1]  <- a0 + chol(P0)%*%rnorm(ne)
  zeta      <- mvtnorm::rmvnorm(T+npad, mean = rep(0,ne), sigma = psi)# Latent variable residuals.
  epsilon   <- rmvnorm(T+npad, mean = rep(0,ny), sigma = theta)  # Measurement errors.
  # Generate latent factor scores
  for (i in 2:(T+npad)){
    etaC[ ,i]   <- eta_intercept + beta %*% etaC[ ,i-1] +
                  G%*%zt[(i+(p-1)*T-1),] + zeta[i, ] + shocks[(i+(p-1)*T),]
  }#End of loop over T
  eta <- t(etaC[,(npad+1):(npad+T)])
  # generate observed series
  y   <- matrix(0, nrow = T, ncol = ny)
  for (i in 1:nrow(y)){
    y[i, ] <- lambda %*% eta[i, ] + epsilon[i, ]
  }
  yall[(1+(p-1)*T):(p*T),] = y
  etaall[(1+(p-1)*T):(p*T),] = eta
} #End of p loop

yall = data.frame(ID = rep(1:np,each=T),Time = rep(1:T,np),
                  V1 = yall[,1], V2 = yall[,2],z1=zt)
etaall = cbind(rep(1:np,each=T),rep(1:T,np),etaall)
colnames(etaall) = c("ID","time",paste0("LV",1:ne))
etaall = data.frame(etaall)


mytheme <- theme_classic() %+replace% 
  theme(axis.title.x = element_blank(), 
        axis.title.y = element_text(face="bold",angle=90))  
p1 <- ggplot(data = etaall, aes(x = time, y = LV1, group = ID, colour = ID)) +
  mytheme +
  labs(title = "Simulated LV", y = "LV 1", x = "Time") + 
  geom_line(size=1) + theme(legend.position="none")

p2 <- ggplot(data = etaall, aes(x = time, y = LV2, group = ID, colour = ID)) +
  mytheme +
  labs(title = "Simulated LV", y = "LV 2", x = "Time") + 
  geom_line(size=1) + theme(legend.position="none")

p3 <- ggplot(data = yall, aes(x = Time, y = V1, group = ID, colour = ID)) +
  mytheme +
  labs(title = "Simulated data", y = "Observed y1", x = "Time") + 
  geom_line(size=1) + theme(legend.position="none")

p4 <- ggplot(data = yall, aes(x = Time, y = V2, group = ID, colour = ID)) +
  mytheme +
  labs(title = "Simulated data", y = "Observed y2", x = "Time") + 
  geom_line(size=1) + theme(legend.position="none")

grid.arrange(p1, p2, p3, p4, ncol = 2)

if (!is.null(fileObs)){
write.table(yall,fileObs,row.names=FALSE,col.names=FALSE,sep=",")}
if (!is.null(fileState)){
write.table(etaall,fileState,row.names=FALSE,col.names=FALSE,sep=",")}

return(list(yall = yall, etaall = etaall))
}

4 Generate data using a VAR(1)-X model

In this tutorial, we use a group-based vector autoregressive model of order 1 (lag 1) with one exogenous input variable, \(z_i(t)\). That is, we are using a VAR(1)-X model, with both process noises and measurement errors, and we assume that all units (e.g., individuals) within the group can be described by the same model with the same parameters.

Further more, we implement a scenario where the the system’s desired target level is at the value of 0, but we incur unusually large shocks at some randomly selected, individual-specific time points such that the system is “perturbed” and pushed away even further from its target level than it normally would be. We show that by applying the optimal input values calculated by the controller*, the system can return to its target level more quickly than if we just let the system run its course and “recover” at its own rate.

The VAR-X model is specified in state-space form as: \[ \text{Measurement model:}\\ V_{1i}(t) = \eta_{1i}(t) + \epsilon_{1i}(t)\\ V_{2i}(t) = \eta_{2i}(t) + \epsilon_{2i}(t)\\ [\epsilon_{1i}(t), \epsilon_{2i}(t)]' \sim N([0,0]',diag(var_{e1},var_{e2}))\\ \text{Dynamic model:}\\ \eta_{1i}(t) = \phi_{11}\eta_{1i}(t-1) + \phi_{12}\eta_{2i}(t-1) + g_1 z_{i}(t-1) + \zeta_{1i}(t)\\ \eta_{2i}(t) = \phi_{21}\eta_{1i}(t-1) + \phi_{22}\eta_{2i}(t-1) + g_2 z_{i}(t-1) + \zeta_{2i}(t)\\ [\zeta_{1i}(t), \zeta_{2i}(t)]' \sim N\left([0,0]', \left(\begin{array}{cc} \psi_{11} & \psi_{12}\\ \psi_{12} & \psi_{22} \end{array}\right) \right) \]

5 Specify VAR(1)-X model in dynr and cook it.

Specify the dynr code to fit the VAR(1)-X model

## Read in data and set up data structure in dynr
etaall = read.table(file1,header=FALSE,sep=",")
thedat = read.table(file0,header=FALSE,sep=",")
colnames(thedat) = c("ID","Time",paste0("V",1:ny),"z1")

## Manually shift the covariates because in dynr, eta_t is matched with u_t in the
#dynamic model
thedat <- plyr::ddply(thedat, .(ID), mutate,
                      zlag1 = dplyr::lag(z1,1,default=0))

## Set up dynr Data
thedat2 <- dynr.data(thedat, id="ID", time="Time", 
                     observed=paste0("V",1:ny),covariates="zlag1")

## VAR(1) recipe specification code starts here

## Prepare the dynamic model
formula=list(
  eta1~ phi11*eta1 + phi12*eta2 + g1*zlag1,
  eta2~ phi21*eta1 + phi22*eta2 + g2*zlag1)

dynamics<-prep.formulaDynamics(formula=formula,
                               startval=c(phi11=.5,phi22=.5,
                                          phi12=-.1,phi21=-.1,
                                          g1=-.5, g2=-.5),
                               isContinuousTime=FALSE) 

## Meausurement model
meas <- prep.measurement(
  values.load=matrix(c(1,0,
                       0,1),ncol=ne,byrow=TRUE), #Starting values for entries in Lambda
  params.load=matrix(c(rep('fixed',ne),
                       rep('fixed',ne)),
                     ncol=ne,byrow=TRUE), #Labels for fixed and freed parameters 
  state.names=c("eta1","eta2"), #Labels for latent variables in eta(t)
  obs.names=paste0('V',1:ny) #Labels for observed variables in y(t)
)

## Initial condition means and covariance matrix. Here we specify them as fixed.
initial <- prep.initial(
  values.inistate=c(0, 0),
  params.inistate=c('fixed','fixed'),
  values.inicov=matrix(c(1,0,
                         0,1),ncol=ne,byrow=TRUE),
  params.inicov=matrix(c('fixed',0,
                         0,'fixed'),ncol=ne))


## Process and measurement noise covariance matrices
mdcov <- prep.noise(
  values.latent=matrix(c(2,.5,
                         .5,1.5),ncol=ne,byrow=TRUE),
  params.latent=matrix(c('psi_11','psi_12',
                         'psi_12','psi_22'),ncol=ne,byrow=TRUE), 
  values.observed=diag(rep(.5,ny),ny), 
  params.observed=diag(paste0('var_e',1:ny),ny)
)

## Put recipes and data together to prepare the full model
model <- dynr.model(dynamics=dynamics, measurement=meas,
                    noise=mdcov, initial=initial, data=thedat2,
                    outfile="VAR.c")
## Cook it!
res <- dynr.cook(model,debug_flag=TRUE,verbose = FALSE)
## Optimization function called.
## Starting Hessian calculation ...
## Finished Hessian calculation.
## Original exit flag:  3 
## Modified exit flag:  3 
## Optimization terminated successfully: ftol_rel or ftol_abs was reached. 
## Original fitted parameters:  0.705573 0.593608 -0.3030289 -0.2029258 -0.7002757 
## -0.4004068 0.6709134 0.2405456 0.3763576 -0.6645845 -0.7833649 
## 
## Transformed fitted parameters:  0.705573 0.593608 -0.3030289 -0.2029258 
## -0.7002757 -0.4004068 1.956023 0.4705127 1.570148 0.5144873 0.4568661 
## 
## Doing end processing
## Successful trial
## Total Time: 1.967817 
## Backend Time: 1.957187
##      phi11      phi22      phi12      phi21         g1         g2     psi_11 
##  0.7055730  0.5936080 -0.3030289 -0.2029258 -0.7002757 -0.4004068  1.9560232 
##     psi_12     psi_22     var_e1     var_e2 
##  0.4705127  1.5701478  0.5144873  0.4568661
## Coefficients:
##         Estimate Std. Error t value  ci.lower  ci.upper            Pr(>|t|)    
## phi11   0.705573   0.006487  108.78  0.692860  0.718286 <0.0000000000000002 ***
## phi22   0.593608   0.008787   67.56  0.576386  0.610829 <0.0000000000000002 ***
## phi12  -0.303029   0.008334  -36.36 -0.319364 -0.286694 <0.0000000000000002 ***
## phi21  -0.202926   0.005263  -38.56 -0.213240 -0.192611 <0.0000000000000002 ***
## g1     -0.700276   0.013440  -52.10 -0.726618 -0.673934 <0.0000000000000002 ***
## g2     -0.400407   0.012059  -33.20 -0.424042 -0.376772 <0.0000000000000002 ***
## psi_11  1.956023   0.052265   37.42  1.853585  2.058461 <0.0000000000000002 ***
## psi_12  0.470513   0.019835   23.72  0.431637  0.509388 <0.0000000000000002 ***
## psi_22  1.570148   0.047813   32.84  1.476435  1.663860 <0.0000000000000002 ***
## var_e1  0.514487   0.035314   14.57  0.445273  0.583701 <0.0000000000000002 ***
## var_e2  0.456866   0.035128   13.01  0.388016  0.525716 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log-likelihood value at convergence = 111046.08
## AIC = 111068.08
## BIC = 111151.86

6 Control theory input computation

Here we proceed as follows:

  1. We simulate another set of data with the true parameters.

  2. Add shocks at particular time points. The shocked latent variables are in X_noControl.

  3. To implement the control theory algorithm, we set the parameters to those estimated with dynr using the previous test data set (primarily \(B\) and \(G\)). and compute the optimal control input.

  4. We need to set the design matrices (\(Q\), \(R\), \(Q_f\)), and control horizon, \(h\).

  5. We use \(E(\eta_i(t)|y_i(1),y_i(2),\ldots,y_i(t-1))\) or \(E(\eta|y_i(1),y_i(2),\ldots,y_i(T_i))\) from dynr to replace X_noControl.

7 Applying LQ control to the new data set.

The basic concept of RHC can be summarized as follows. At the current time \(t\), the optimal control is computed, over a finite fixed horizon, \([t,t+h]\), that is, from time points \(t\), \(t+1, \ldots, t+h\). Of the control values over this time window, only the one associated with time t is adopted and administered to control the state values of the system at time \(t + 1\). The procedure is then repeated sequentially for the next time, namely, for time t+1, over the horizon \([t+1; t+1+h]\), to yield a declining (receding) possible horizon over which future control values may be administered before the end of the time series.

You may comment out the two for loops over Rmult and h_cond to repeat the control theory computations for a range of control horizon values (h = 4, 10) and different cost values for administering the control input (R = 0.001, 0.1, 1)

Rmult = 0.001; h_cond = 4
#for (Rmult in c(0.001, .1, 1)){
#  for (h_cond in c(4,10)){

# Define Q and R
# Q =  is a matrix of cost efficients that expresses the penalty to be imposed on 
# state deviations from the desired level.
# The maximum values depend on the values of X. Could set to identity matrix or
# t(lambda)%*%lambda, and play with the values of R.

# R = is a cost matrix for penalizing against the administration of input. 
# Usually start with something small, like .1

Q = t(lambda)%*%lambda   #ne x ne weight matrix for state deviations
Q_f = t(lambda)%*%lambda #ne x ne weight matrix for state deviations
                         #at the terminal time point of each control horizon
R = Rmult*diag(1,r) #Weight matrix of size r x r representing cost of 
                    #administering control input
I_q = diag(1,ne) #Identity matrix of size ne
h = h_cond #control horizon
eta_target = matrix(0,T,ne)#Same set of target (reference) levels of states
                           #constrained to be the same across individuals.
                           #Could be made person-specific.

X_NoControl = matrix(NA, ncol = ne, nrow = np*T) #Place holder for latent variables
X_dev  = matrix(NA, ncol = 5, nrow = np*T) #Place holder to keep track of 
                                           #discrepancies between state variable                                               #values and target levels
thedat4 = thedat3
for (i in 1:np){
  temp.dat = thedat3[(1+(i-1)*T):(i*T),]
  #Set up dynr data and model
  thedat.original = dynr.data(temp.dat, id="ID", time="Time", 
                              observed=paste0("V",1:ny),covariates="zlag1")
  model.Nocontrol <- dynr.model(dynamics=dynamics, measurement=meas,
                                noise=mdcov, initial=initial, data=thedat.original,
                                outfile="temp.c")
  coef(model.Nocontrol) <- coef(res)
  res.Nocontrol = dynr.cook(model.Nocontrol,optimization_flag = FALSE,
                            hessian_flag = FALSE,
                            verbose = FALSE, debug_flag=TRUE)
  X_Originali = res.Nocontrol@eta_filtered
  X_NoControl[(1+(i-1)*T):(i*T),] = as.matrix(t(X_Originali),byrow=TRUE,
                                              ncol=ne)
  temp.dat$oldzlag1 = temp.dat$zlag1
  
  L <-list(matrix(0,ncol=ne,nrow=ne),h)
  g <- list(matrix(0,ncol=1,nrow=ne),h)
  zstar <- matrix(0,nrow=r,ncol=1)
  
  L[[h]] <- Q_f  # 7a
  g[[h]] <- -Q_f %*% eta_target[1,] #7b
  for (t in 1:(T-1)){
    #  thedat.test = dynr.data(temp.dat[1:(min(T,t+h)+(i-1)*T),], id="ID", time="Time", 
    #                              observed=paste0("V",1:ny),covariates="oldzlag1")
    #  model.test <- dynr.model(dynamics=dynamics, measurement=meas,
    #                                noise=mdcov, initial=initial, data=thedat.test,
    #                                outfile="temp0.c")
    #  coef(model.test) <- coef(res)
    #  res.test = dynr.cook(model.test,optimization_flag = FALSE,
    #                            hessian_flag = FALSE,
    #                            verbose = FALSE, debug_flag=TRUE)
    #  X_test = res.test@eta_filtered #Set XOriginalitPlus to elements in X_test if done recursively
    for (tau in (h-1):1) {
      inputMultiplier <- solve( I_q + (L[[tau+1]] %*% G2 %*% solve(R)%*% t(G2))) #From Kwan & Han
      L[[tau]] <- t(B) %*% inputMultiplier%*% L[[tau+1]] %*%B + Q #7e from Molenaar; Eq 26, Kwan & Han
      g[[tau]] <- t(B) %*% inputMultiplier%*% g[[tau+1]] - Q %*% eta_target[min(t+tau,T),] #7f
      if (tau==1){
        X_OriginalitPlus = matrix(X_Originali[,thedat.original$time ==t+tau],ncol=1)
        zstar <- -solve(R) %*% t(G2) %*% inputMultiplier %*% (L[[tau+1]] %*% B %*% 
                                                                X_OriginalitPlus + g[[tau+1]]) #7d; see K&H, 5.67, top, p. 231        
      }}
    temp.dat$z1[t] = zstar 
    temp.dat$zlag1[t+1] = zstar #usually this should be [t], set to [t+1] because of dynr setting for cov in the dy eqs.
  } #T-loop
  thedat4[(1+(i-1)*T):(i*T),] = temp.dat[,!colnames(temp.dat) %in% "oldzlag1"]
}

#One more pass to get everyone's smoothed state estimates with zstar
thedat.Controlled = dynr.data(thedat4, id="ID", time="Time", 
                              observed=paste0("V",1:ny),covariates="zlag1")
model.Controlled <- dynr.model(dynamics=dynamics, measurement=meas,
                               noise=mdcov, initial=initial, data=thedat.Controlled,
                               outfile="temp.c")
coef(model.Controlled) <- coef(res)
res.Controlled <- dynr.cook(model.Controlled,optimization_flag = FALSE,
                            hessian_flag = FALSE,
                            verbose = FALSE, debug_flag=TRUE)
X_Controlleds <- as.matrix(t(res.Controlled@eta_filtered),
                           ncol=ne,byrow=TRUE) #Extract filtered estimates

#Extract zstar
zControlled = as.matrix(thedat4$z1,col=1)

#Below we are going to regenerate simulated data with the same random number seed, but now using zstar.
set.seed(12345)
controlData3 = GenerateData(T,npad,np,ne,ny,r,psi,lambda,theta,beta,eta_intercept,
                            G,a0,P0,zControlled,
                            fileObs="Controlled.obs", fileState="Controlled",shocks)

8 Evaluate output

Here, we inspect some summary statistics and plots to evaluate the results of applying the controller vs. not.

##   True_NoShock     True_Shock     True_zstar   Filtered_old Filtered_zstar 
##       1.970193       3.014143       2.480108       2.755830       2.580061

Notice that with the control input (under “True_zstar”), the system shows smaller root mean squared deviation from the target level than when no control input is applied (see “True_Shock”), and the root mean squared deviation value is just slightly higher than if no shocks are applied to the system. SImilarly, the filtered estimates with control input (“Filtered_zstar”) also have smaller root mean squared deviation value than the filtered estimates with no control (“Filtered_old”).

Verify from the plots over time that when the control input values are applied, the individual trajectories within the system all return to their target baseline (\(=\) 0) more efficaciously than if no control input is fed into the system.

np_test = c(1,2,3)
for (id in np_test){
for (toPlot in 1:ne){
  #pdf(paste0("ControlChapterh",h,"R",R,"ID",id,"X",toPlot,".pdf"))
  par(cex.axis=1.2,cex.lab=1.3,cex.main=1.4,mar=c(2.1,4.1,2.1,4.1))
  plot(1:T,X_Controlled[1:T +(id-1)*T,toPlot], 
       xlim=c(0,T),
       ylim = c(min(X_Shocks[,toPlot],na.rm=T)-2,                                               max(X_Shocks[,toPlot],na.rm=T)+.1),
     ylab=bquote(eta[.(toPlot)]),
     xlab="Time",
     main=paste0("ID = ",id,', R = ',R,", h = ",h, ' with and w/o control'),
     type="n")
lines(1:T,X_Controlled[1:T +(id-1)*T,toPlot],col="blue")
points(1:T,X_Controlled[1:T +(id-1)*T,toPlot],col="blue",pch="C")
lines(1:T,X_Shocks[1:T +(id-1)*T,toPlot],col="red")
points(1:T,X_Shocks[1:T +(id-1)*T,toPlot],col="red",pch=5)
abline(h=eta_target[id,toPlot], col = 'black', lty =2 )
seqq = seq(min(X_Shocks[,toPlot]-2,na.rm=T),            max(X_Shocks[,toPlot],na.rm=T),1)
lines(rep(shockPoints[id,1],length(seqq)),seqq,col="gray",lwd=3)
lines(rep(shockPoints[id,2],length(seqq)),seqq,col="gray",lwd=3)
legend("top", legend=c(as.expression(bquote(eta[.(toPlot)~"controlled"])),
                       as.expression(bquote(eta[.(toPlot)~"no control"])),
                       "Shock points",
                       as.expression(bquote("Target level,"~eta^r)),
                       "Control input z*"),
         col=c("blue", "red","gray","black","darkgreen"), lty=c(1,1,1,2,3), 
         cex=1,pch=c("C","NC",NA,NA,"z*"),bty="n",lwd=c(1,1,3,1,1))
par(new=TRUE,mar=c(2.1,4.1,2.1,4.1))
plot(1:T,thedat4$zlag1[thedat4$ID==id],type="n",
    xlab=NA,ylab=NA,xlim=c(0,T),
    ylim=c(-6, 14),axes=FALSE)
lines(1:T,thedat4$zlag1[thedat4$ID==id],col="darkgreen",lty=3)
points(1:T,thedat4$zlag1[thedat4$ID==id],pch="z",col="darkgreen")
axis(side = 4)
mtext(text='Control input (z*)',line=2.5,
      side = 4, at = 4.5, 
      outer=FALSE,cex=1.3)
#dev.off()
}}