OPT_DES
Optimal controller design based upon the methodology of Akaike and Nakagaw (1972). The routine allows for multiple channels for both the controlled and manipulated variables. The gain matrix is computed after automatically selecting the best autoregressive model using minimum multivariate final prediction error (MFPE).
Required Arguments
MAXLAG — Maximum number of autoregressive parameters requested. (Input)
X — NOBS by NCHANX matrix containing the multi‑channel time series for the manipulated variables, also called control system input variables. Each row corresponds to a different observation in the series, and each column of X corresponds to a univariate time series for one of the controlled channels. (Input)
Y — NOBS by NCHANY matrix containing the multi‑channel time series for the controlled variables, also called control system output variables. Each row corresponds to a different observations in the series, and each column of Y corresponds to a univariate time series for one of the manipulated channels. (Input)
NPAR — Number of autoregressive parameter matrices in minimum FPE model. (Output)
GAIN — NCHANX by NPAR gain matrix. (Output)
Optional Arguments
Q — NCHANY by NCHANY non‑negative definite symmetric weighting matrix for the quadratic optimization criterion. (Input)
Default: If AVAR is non‑singular, Q = Inverse of AVAR. If AVAR is singular, then
Q = diag(1/AVAR1,1, 1/AVAR2,2, …, 1/AVARNCHANY,NCHANY).
R — A NCHANX by NCHANX positive definite symmetric weighting matrix for the quadratic optimization criterion. (Input)
Default: R = diag(1/σ2X1, 1/σ2X2, …, 1/σ2XNCHANX)
IPRINT — Printing option. (Input)
|
0 |
No printing |
|
1 |
Prints final results only |
|
2 |
Prints intermediate and final results |
Default: IPRINT = 0.
NCHANX — Number of time series for manipulated variables. (Input)
Default: NCHANX = size(X,2)
NCHANY — Number of time series for controlled variables. (Input)
Default: NCHANY = size(Y,2)
NOBS — Number of observations in each time series. (Input)
Default: NOBS = size(X,1)
NSTAGE — Number of stages used to compute the gain matrix, GAIN. (Input)
Default: NSTAGE = 20
TRANSY —NPAR by NCHANY by NCHANY transition matrix. (Output)
GAMMAX — NPAR by NCHANY by NCHANX gamma matrix. (Output)
CCV — NCHAN by NCHAN by (MAXLAG+1) matrix containing the sample autocovariances of the NCHAN time series variables, where NCHAN = NCHANX+NCHANY. For the i‑th time series variable, the first element, CCV(i,i,1) is the sample variance for the i‑th series and the remaining elements, CCV(i,i,2) … CCV(i,i,MAXLAG+1), contain the autocovariances of the i‑th series for lags 1 thru MAXLAG. Elements CCV(i,j,2) … CCV(i,j,MAXLAG +1) contain the autocovariances between the i‑th and j‑th series for lags 1 through MAXLAG. (Output)
AVAR — NCHANY by NCHANY matrix containing estimates of the noise variances, autocovariances and cross-covariances for the NCHANY time series. (Output)
FPE — Final Prediction Error for fitted model. (Output)
AIC — Akaike’s Information Criterion for fitted model. (Output)
FORTRAN 90 Interface
Generic: CALL OPT_DES (MAXLAG, X, Y, NPAR, GAIN [, …])
Specific: The specific interface names are S_OPT_DES and D_OPT_DES.
Description
The routine OPT_DES is based upon the FPEC and OPTDES program published in the TIMSAC –71 Library described by Akaike, H. and Nakagawa, T (1972). Estimates of the autoregressive parameters for the model with minimum multivariate final prediction error (MFPE) are calculated using the methodology described in Akaike, H., et. al (1979). Their methodology produces an estimate of the gain matrix, G, in a feedback control system with the following relationship:
where
Xt is a vector containing the value of the NCHANX control variables at time = t,
Zt = the state vector at time = t, and
G is the NCHANX by NPAR gain matrix (GAIN).
The gain matrix is estimated by minimizing the quadratic criterion:
where E is the expectation operator, and I = NSTAGE.
Akaike and Nikagawa (1972) describe a process for obtaining the optimum gain matrix by adjusting the quadratic criterion matrices, Q and R. Initially they recommend setting Q and R to the default values described above for Q and R. They recommend that the behavior of the controller using the gain matrix obtained from these defaults be examined by simulating the controller.
From these simulations, the variances of the control input variables (X) should be examined to identify channels whose variances are too large or too small. If they are too large or small, their corresponding diagonal elements in R and Q should be decreased or increased.
The multivariate final prediction error for a multivariate autoregressive model with lag p is defined as:
where
is the determinant of the estimated NCHANY × NCHANY matrix of covariances of the white noise in the multivariate series, N =NOBS and m =NCHANY.
The model selected and parameter estimates vary depending upon the value of MAXLAG. Akaike and Nakagawa (1972) provide a rule of thumb for MAXLAG: start with values between 
and 
, and keep MAXLAG below 
.
Similar to the univariate case, MAXLAG must never exceed 
.
The numerical accuracy decreases as MAXLAG increases. In this case, it is possible for the estimated MFPE to become negative. If this happens try using double precision.
NSTAGE, the number of stages used to compute the gain matrix, should be selected to provide an accurate estimate of the gain matrix. Essentially, the gain matrix is the matrix for control of the system at time = NSTAGE. As NSTAGE is increased, the gain matrix converges to a constant solution. That is, as NSTAGE is increased,
It is tempting to set NSTAGE to a large value. However, this will increase execution time and subject the final estimates to rounding errors.
Example
The following example uses the Gas Furnace Data (Box and Jenkins 1976, pages 532‑533). In this example, the controller has one variable, percent CO2 , that is controlled by a single manipulated variable, input gas rate in cubic feet/minute. These multi‑channel series consist of NOBS=296 observations.
USE GDATA_INT
USE WRRRN_INT
USE OPT_DES_INT
IMPLICIT NONE
INTEGER, PARAMETER :: MAXLAG=10, NOBS=296
INTEGER NCOL, NROW, NPAR
REAL(KIND(1E0)) RDATA(NOBS,2), X(NOBS,1), Y(NOBS,1), GAIN(1,MAXLAG), &
Q(1,1), R(1,1)
EQUIVALENCE (X, RDATA(1,1)), (Y, RDATA(1,2))
Q(1,1) = 0.16E0
R(1,1) = 0.001E0
GAIN=0.0E0
CALL GDATA (7, RDATA, NROW, NCOL)
CALL OPT_DES(MAXLAG, X, Y, NPAR, GAIN, Q=Q, R=R)
CALL WRRRN("GAIN", GAIN(:,:NPAR))
END
Output:
GAIN
1 2 3 4
1.418 1.851 1.894 1.502
