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