This document summarizes the contents of N. L. Ricker's Tennessee Eastman
Challenge Problem archive.
TEPROB.F Fortran code provided by Tennessee Eastman (ASCII FORMAT)
TECOMMON.INC An "include" file needed to compile TEPROB.F (ASCII FORMAT)
TE_MEX.F Fortran code needed to generate the .MEX interface
to Matlab. (ASCII FORMAT)
TE_MEX_TC.FV4 was written for Matlab 4.1 on a MacIntosh. It
has not been tested on other machines. It is similar to TE_MEX.F, but includes
PI control of the reactor temperature and separator temperature. There are
also modifications to the output vector. See comments in the code for more
information. The state vector for this version has 52 elements. The last 2 are the
integrated errors for the 2 PI controllers. You must initialize these properly
for the system to be at steady-state.
The following 3 files are in tab-delimited text (ASCII) format.
It should be possible to read them into a spreadsheet program or text editor.
They contain both text (row labels and units) and numerical data, so you may
need to re-format before reading into Matlab, etc. The files give operating
conditions at 6 different steady-state operating conditions described in
the paper "Optimal Steady-state Operation of the Tennessee Eastman
Challenge Process", N. L. Ricker, Computers & Chemical Engineering, Vol. 19,
No. 9, pp. 949-959(1995). NOTE: in revision, Table 2 was eliminated. It
gives the first 38 state variables at the base case and the 6 optimal
operating modes. Tables 3 and 4 were renumbered 2 and 3, respectively.
To generate the complete state vector (50) variables, you must append
the 12 manipulated variables (from Table 4) to the end of the 38 states
in Table 2. To make sure that it was done correctly, verify that it is a
steady-state condition by calling TEFUNC to calculate the derivatives
for the given states and manipulated variables. All derivatives should
be less than 1.E-3.
TABLE2.TXT Table 2 -- state variables
TABLE3.TXT Table 3 -- output variables
TABLE4.TXT Table 4 -- manipulated variables
The following 3 directories contain Matlab and Fortran code described in
the papers "Nonlinear Model Predictive Control of the Tennessee Eastman
Challenge Process", Computers & Chemical Engineering, Vol. 19, No. 9,
pp. 961-981(1995), and "Nonlinear Modeling and State Estimation for the
Tennessee EastmanChallenge Process", ibid, pp. 983-1005(1995).
TEest3 Directory containing code related to the above papers. See
the file LARRY$DUA1:[ricker.TE_send.TEest3]TEest3.doc
for more information on these files.
TEest6 Directory containing code related to a newer version of the TEest3
model. See the file LARRY$DUA1:[ricker.TE_send.TEest6]example.doc
for more information.
nl_mpc Code to implement nonlinear MPC. See CONTENTS.TXT (in directory)
for more information.
The following directory contains Matlab and Fortran code described in the
paper "MPC of a continuous, nonlinear, two-phase reactor", N. L. Ricker,
J. Process Control, vol. 3, 109-123(1993).
TE4
See the file contents.txt for more information on the contents of the TE4
directory.
IDV1 Directory containing data files from a simulation of the above
decentralized strategy with disturbance 1 of Downs and Vogel.
The disturbance occurs at time = 1 hour and lasts for the
remainder of the 50-hour simulation. Variables are sampled
every 10 minutes. There are 4 data files, organized as follows:
y.dat ASCII (TEXT) file containing the measured outputs. Each row
represents a time instant, starting at time = 0. Each column
is a measured variable. The first 41 are the XMEAS variables
from the TE model (measured outputs). Columns are delimited by
TAB characters to make it easier to load the files into plotting
and analysis packages (such as Matlab). The remaining columns
are as follows:
(42) is for cost [cents/kmol product].
(43) is production rate of G [kmol G generated/h]
(44) is production rate of H [kmol H generated/h]
(45) is production rate of F [kmol F generated/h]
(46) is partial pressure of A in reactor [kPa]
(47) is partial pressure of C in reactor [kPa]
(48) is partial pressure of D in reactor [kPa]
(49) is partial pressure of E in reactor [kPa]
(50) is true (delay free) mole % G in product
(51) is true (delay free) mole % H in product
These "extra" outputs are not used in the control strategy.
They may be useful in interpreting the results, however.
u.dat As for y, but containing the 12 manipulated variable signals
(XMV) of the TE model. All values are in the range 0-100 %.
r.dat As for y, but containing 36 variables used in the decentralized
control strategy. These include setpoints for the various loops
as well as other signals. They are probably not of general
interest, but just in case, here's a description.
(1) Setpoint for A feed flow (stream 1), kscmh
(2) Setpoint for D feed flow (stream 2), kg/h
(3) Setpoint for E feed flow (stream 3), kg/h
(4) Setpoint for C+A feed flow (stream 4), kscmh
(5) Setpoint for purge rate (stream 9), kscmh
(6) Setpoint for sep. underflow (stream 10), m^3/h
(7) Setpoint for product rate (stream 11), m^3/h
(8) Setpoint for reactor pressure, kPa
(9) Setpoint for reactor level, %
(10) Setpoint for separator level, %
(11) Setpoint for stripper level, %
(12) No longer used
(13) Production rate target (stream 11), m^3/h, supplied
by the operator. Note that this may be overridden under
some conditions. Actual setpoint for stream 11 is in
(7), above.
(14) Target for mole % G in product (stream 11), as supplied
by the operator. May be modified (rate-of-change constraint).
The setpoint used in the feedback loop is (34), below.
(15) No longer used
(16) Setpoint for %A/(%A + %C) in reactor feed (stream 6), %
(17) Setpoint for %A + %C in reactor feed (stream 6), %
(18) No longer used
(19) Maximum reactor pressure -- setpoint for pressure override,
Loop 18 as described in the paper (kPa).
(20) Minimum value for separator coolant valve (reactor level
override control), %. See discussion of Loop 19 in the paper.
(21) Maximum value for separator coolant valve (reactor level
override control), %. See discussion of Loop 19 in the paper.
(22) Ratio setpoint defining variable (1) = (22) * (32).
(23) Ratio setpoint defining variable (2) = (23) * (32).
(24) Ratio setpoint defining variable (3) = (24) * (32).
(25) Ratio setpoint defining variable (4) = (25) * (32).
(26) Ratio setpoint defining variable (5) = (26) * [(32) - (33)].
(27) Ratio setpoint defining variable (6) = (27) * [(32) - (33)].
(28) Ratio setpoint defining variable (7) = (28) * [(32) - (33)].
(29) No longer used
(30) Output of reactor pressure override loop. Used to decrease
production rate when reactor pressure is too high.
(31) Production rate index. Nominal value = 100*(var 13)/23. May
be modified to limit rate-of-change.
(32) Production rate index after adjusting for overrides, contraint
on rate of change, and feedback from production rate loop.
Used to determine setpoints 1-7. See formulas for variables
(22)-(28), above.
(33) Feedback adjustment to production rate index (Loop 8 of paper).
(34) Current setpoint used in feedback control of mol % G in
stream 11.
(35) Eadj value used in eqs. 5 and 6 of the paper.
(36) Reactor temperature setpoint, degrees C.
t.dat The first column gives the time instants corresponding to
each row in y, u and r. The other 2 columns in t.dat are
not of general interest and can be ignored.
IDVx where x is a number between 1 and 20 is as for IDV1.dir, but for
disturbance #x of Downs and Vogel.
USING THE MEX FILE IN MATLAB:
For general information on S-functions, see the SIMULINK manual.
In addition to the standard input arguments (time, states, inputs, flag),
TE_MEX requires 4 parameters, as follows:
p1 idv(20) -- a vector of length 20 containing the desired values of
the disturbance flags (see paper by Downs and Vogel).
The following are all used by the TE routines to save data between
Matlab calls. They are initialized by a call to this function with
time<0 and flag=0. You must create Matlab variables of the
appropriate size (a vector of zeros is convenient),
and pass them to te_mex as input arguments. Once you have
initialized these variables, DO NOT modify them during
the simulated plant operation.
p2 Vector of length 153
p3 Vector of length 586
p4 Vector of length 139
Here are some example calls from Matlab:
t=-1;
x=[];
u=[];
flag=0;
idv=zeros(20,1);
p2=zeros(153,1);
p3=zeros(586,1);
p4=zeros(139,1);
[sys,x0]=te_mex(t,x,u,flag,idv,p2,p3,p4);
u=[specify a vector, length 12];
flag=1;
t=0;
dxdt=te_mex(t,x0,u,flag,idv,p2,p3,p4);
flag=3;
y=te_mex(t,x0,u,flag,idv,p2,p3,p4);
The first call to te_mex (with t=-1) initializes the storage vectors, p2 to p4,
and returns the (default) initial state vector in x0. This is a good test
to verify that the code is working (especially for those of you who have
compiled and installed it yourselves). You can compare x0 to the values
listed for YY(1) to YY(50) in the TEINIT subroutine supplied
by Downs and Vogel. Doing this, for example, helped me to figure out how
to port the code from a Mac, where it was developed, to a VAX. Our installation
requires the G_Float option (see Matlab/VAX manual), and the x0 vector was
incorrect unless this was set properly for the compilation
(and the fmexg .com file was used for linking).
The second call to te_mex (with flag=1) calculates the time-derivatives of
the states at the initial conditions (x0) and for a set of inputs specified by
the "u" vector. The derivatives are returned in the "dxdt" variable.
The third call (with flag=3) returns the values of the outputs (all 41 of them -
see paper by Downs and Vogel) for the given x0 and u. As explained in the
paper, when t>0, measurement noise is added to the outputs automatically.
Also, the composition analyzers are only sampled at certain times. Thus,
proper specification of the time variable is important.
The typical use of TE_MEX.MEX is to simulate operation of the plant for
a given set of conditions.
For example, the following sequence of Matlab statements
shows how to reproduce the simulations of Figure 3 of Downs and Vogel.
The SIMULINK function, RK45, is used to integrate the differential
equations for a period of 0.3 hours. See the SIMULINK manual for details
on this function.
% Initialize variables:
idv=zeros(20,1);
p2=zeros(153,1);
p3=zeros(586,1);
p4=zeros(139,1);
[sys,x0]=te_mex(-1,[],[],0,idv,p2,p3,p4);
u0=[63.53, 53.98, 24.644, 61.302, 22.21, ...
40.064, 38.1, 46.534, 47.446, 41.106, ...
18.114, 50]; % Base case inputs (XMV variables).
u0(10)=38; % Specifies a step in XMV(10)
% Set up for integration:
tvec=[0; 0.3]; % Starting and ending time.
ut=[tvec [u0;u0] ]; % Specifies constant inputs.
options=[.001, .00001, .001, 0, 0, 0, 0];
% Integrate over entire 0.3 hour time period:
[tt,xt,yt]=rk45('te_mex',tvec,x0,options,ut,idv,p2,p3,p4);
% Plot results:
clf
subplot(221)
plot(tt,yt(:,6)),title('Reactor Feed'),xlabel('Time (hours)')
subplot(222)
plot(tt,yt(:,7)),title('Reac Pressure'),xlabel('Time (hours)')
subplot(223)
plot(tt,yt(:,8)),title('Reac Level'),xlabel('Time (hours)')
subplot(224)
plot(tt,yt(:,9)),title('Reac Temp'),xlabel('Time (hours)')
pause
clf
subplot(221)
plot(tt,yt(:,11)),title('Prod Sep Temp'),xlabel('Time (hours)')
subplot(222)
plot(tt,yt(:,12)),title('Prod Sep Level'),xlabel('Time (hours)')
subplot(223)
plot(tt,yt(:,10)),title('Purge Rate'),xlabel('Time (hours)')
subplot(224)
plot(tt,yt(:,18)),title('Stripper Temp'),xlabel('Time (hours)')
I found that a specification of a maximum step size of about 0.001
was necessary to get reliable accuracy from the RK45 function
(see use of "options" variable, above). Other settings
may give better performance. I didn't experiment with it much.
Elapsed time for the integration was about 1 minute on a Mac IIfx.
The following additional plot commands show how the concentrations
from the analyzers vary with time. You can observe the discrete
sampling behavior, e.g., the product concentrations vary every 0.25
hours, and the others vary every 0.1 hour.
clf
subplot(221)
plot(tt,yt(:,23)),title('A in Feed'),xlabel('Time (hours)')
subplot(222)
plot(tt,yt(:,30)),title('B in Purge'),xlabel('Time (hours)')
subplot(223)
plot(tt,yt(:,40)),title('G in Product'),xlabel('Time (hours)')
subplot(224)
plot(tt,yt(:,41)),title('H in Product'),xlabel('Time (hours)')
Many other uses of TE_MEX.MEX should be obvious to those familiar with
the analysis tools in SIMULINK and the various Matlab toolboxes.
USING THE SOURCE CODE TO COMPILE/LINK YOUR OWN .MEX FILE:
The linking procedure is machine-dependent.
See your Matlab manual. I have verified that the code
works on a Macintosh and a VAX. I would expect it to work on other
machines, but one can never be sure. You're on your own if you
attempt it (or modify the source code).
Some hints follow:
On the Macintosh, I used MPW 3.2 with Language Systems Fortran 3.0. I
compiled TEPROB.F and TE_MEX.F using the following options:
-i '{MexDir}' -mc68020 -mc68881
The symbol {MexDir} was set as follows:
Set MexDir "{Boot}MATLAB:MEX:"
You may need to modify this, depending on your directory structure.
The Link step is tricky. I used the following command:
Link -o "te_mex.mex" -d -f -ad 4 -rt 'MEX0'=0 -sg MEX -m mexmain ¶
-t MEX0 -c MATL -srt -br on ¶
"{MexDir}MPW:cmex.c.o" ¶
"{MexDir}MPW:LS:fmex.c.o" ¶
'te_mex.f.o' ¶
'teprob.f.o' ¶
"{FLibs}FORTRANLib.o" ¶
"{FLibs}IntrinsicLibFPU.o" ¶
"{FLibs}Interface.o" ¶
"{FLibs}Runtime.o"
with the symbol {FLibs} set to "{Boot}MPW:Libraries:Libraries:FLibraries:".
This assumes that you have the compiled object files te_mex.f.o, and
teprob.f.o, stored in the default directory, and it creates te_mex.mex in
the same directory. Modify as appropriate. The key to making it work
is the "-br on" option, which allows linking of large code segments. It
is apparently unavailable in MPW prior to version 3.2.
On a VAX (VMS), rename the .F files to .FOR. Then compile as follows:
$ for/extend/g_float teprob
$ for/extend/g_float te_mex
$ fmexg te_mex,teprob pv,dvec,const,teproc,wlk,randsd
(The fmexg.com file should be available in the Matlab directory on
your VAX system. It links the object code created in the two
compilation steps to generate a file called te_mex.mexg)
N. L. Ricker, Professor
University of Washington
Chemical Engineering
Box 351750
Seattle, WA 98195-1750
E-mail: ricker@cheme.washington.edu
fax: (206) 543-3778
phone: (206) 543-2250 or -8786