By Peter Nachtwey^{1}

Since most of us don’t have our own heat exchanger to test on, I thought it would be good to show how one can be simulated. This would allow us to test the equations that calculate the PI tuning values over the range from a conservative to aggressive response. In addition, the sample period can be changed to see the affect that the sample period has on the response.

**Building the Discrete Time Model
**Learning how to model a first order plus dead time (FOPDT) system is the easiest place to start. This model simply uses the last process variable (PV) and the new CO (control output) to estimate what the new value of PV will be in the next time period.

The form of the discrete time model is:

PV(n+1) = A·PV(n) + B·Kp·CO(n–Өp) + C

The FOPDT model is basically a simple single pole filter similar to an RC circuit in electronics, except here:

1. the system has gain, Kp

2. the input at time n affects the output at time n+1

3. the system has an initial steady state, C, that is not 0

4. the system has dead time, Өp

The first thing that must be done is to build the basic model by specifying the A and B constants. Note that A and B are constants for a first order system and are arrays for higher order systems. So in this study:

A = exp(–T/Tp)

B = (1 – A)

T is the sample period and *T*p is the plant time constant. A and B are coefficients in the range from 0 to 1, and coefficient A tends to be close to 1.

Calculating coefficient C in the above model is the tricky part. The heat exchanger steady state is not 0 °C. If coefficient C is left out the above model equation or it is set to zero, than PV will approach 0. This does not match the observed behavior of the heat exchanger. Therefore, a value for C must be calculated that will provide the correct value of PV when CO is set to 0.

This is done as follows. From the heat exchanger data plot, we can see that the heat exchanger temperature (PV) is 140 °C degrees when the controller output (CO) is 39%. If these values are plugged into the formula y=mx+b, we can solve for the steady state temperature, PVss:

140 = 39Kp + PVss

And this lets us calculate coefficient C:

C = (1 – A) PVss

The data below is from the Validating Our Heat Exchanger FOPDT Model tutorial. The heat exchanger is modeled as a reverse acting FOPDT system with Laplace domain transfer function:

The model parameter values for Kp, *T*p and Өp where provided in the tutorial as:

• Plant gain, Kp = –0.533 °C/%

• Time constant, *T*p = 1.3 min

• Dead time, Өp = 0.8 min

Also,

• Sample time, T = 1/60 min

• Dead time in samples = nӨp = Өp/T = 48

Steady state PV can now be computed assuming the system is linear:

PVss = 140 – 39Kp = 160.8 °C

Next generate the discrete time transition matrix. In this simple case it is a single element:

A = exp(–T/Tp) = 0.987

Also generate the discrete time input coupling matrix, which again is a single element. Notice that the plant gain is negative:

B Kp = (1 – A) Kp

= (1 – 0.987)(–0.533)

= –0.00679

Using the above information, coefficient C becomes:

C = (1 – A) PVss = 2.048

Substitute our values for A, B and C into the discrete time model form above to calculate temperature (PV) as a function of controller output (CO) delayed by dead time, or:

• PV(n+1) = 0.987 PV(n) – 0.00679 CO(n–48) + 2.048

**Open Loop Dynamics
**The easiest way to test this discrete time model is to execute an open loop step test that is similar to the example given in a previous tutorial (see Validating Our Heat Exchanger FOPDT Model).

Initialize:

temperature as: PV(0) = 140

controller output as: CO = 39 for n ≤ 48

Implement controller output step at time t = 25.5 min, or n = 1530:

CO = 39 for n < 1530

CO = 42 for n ≥ 1531

Simulate for 60 min, or: n = 0,1,…,3600. Display results in the plot below after 20 min (click for large view):

**Controller Tuning**

The formulas for calculating the PI tuning values can now be tested using the discrete time process model that was generated above.

The tuning values are for the PI controller that is expressed in the continuous time or analog form as:

The PI controller tuning values are calculated using the equations provided in previous tutorials (see PI Control of the Heat Exchanger). Following that procedure, the desired closed loop time constant, *T*c, must be chosen to determine how quickly or how aggressive the desired response should be.

For an aggressive response, *T*c should be the larger of 0.1*T*p or 0.8Өp. As listed above, *T*p = 1.3 min and Өp = 0.8 min for the heat exchanger. Using these model parameters, then:

Tc = max(0.1*T*p, 0.8Өp) = 0.64

A moderate *T*c of the larger of 1.0*T*p or 8.0Өp will produce a response with no overshoot. We will use this rule in the closed loop simulation below:

Tc = max(1.0*T*p, 8.0Өp) = 6.4

Using this moderate *T*c, we can now compute the controller gain, Kc, and integrator term,*T*i:

**Closed Loop PI Control**

The function logic that follows implements the PI controller and the discrete time FOPDT model’s response. The controller uses the PV to calculate a CO that is used later by the plant model. How much later is determined by the dead time. The model then responds by changing the PV. This new PV is then used in the next iteration.

– compute controller error:

ERR = SP – PV

– compute COp, the proportional contribution to CO:

COp = Kc·ERR

– compute COi, the integral contribution to CO:

COi = COi + (Kc·T·ERR)/*T*i

– compute the total controller output:

CO = COp + COi

– compute the process response at the next step n using the time-delayed CO:

PV(n+1) = A·PV(n) + B·Kp·CO(n–Өp) + C

Aside: The integral contribution to CO can include the integrator limiting technique talked about in the Integral (Reset) Windup and Jacketing Logic post with the expression: |

We initialize the variables at time t=0 by assuming that there is no error and the system is in a quiescent state. Thus, all of the control output, CO, is due only to the integrator term, COi:

PV(0) = 140; CO(0) = 39; COi(0) = CO(0)

The plot below (click for large view) shows the response of the closed loop system over a simulation time of 60 minutes. The set point is stepped from 140 to 138.4 at time t = 25.5 minutes.

The function logic above calculates the PV and CO for each time period in the 60 minute simulation. Note that there are three CO type variables. CO is the current control value, COi is the current integrator contribution to the CO, and CO(n-Өp) is the delayed CO that simulates the dead time.

When programming the functions, additional logic is required keep from using indexes less than 0.

Notice that to prevent overshoot, the rate at which the controller’s integrator output, COi, responds or winds up must be limited so the total control output does not exceed the final steady state control output, and final steady state control output is not reached before the system starts to respond. The way I think about it is that the function is calculating the value for Kc so this goal can be achieved.

Also notice that with moderate tuning, the closed loop response is much slower than the open loop response.

__________

**About the Author**

Peter Nachtwey is the president and head engineer of Delta Computer systems which is a manufacturer of motion control and color sensor products. Peter Nachtwey has more than 20 years of experience developing hydraulic, pneumatic, electrical motion control and vision systems and products for industrial applications.

Peter Nachtwey

President

Delta Computer Systems, Inc.

11719 NE 95th St Suite D

Vancouver, WA 98682-2444

(360) 255-8688

Email: peter@deltamotion.com

Web: http://www.deltamotion.com

**Addendum**

Below is a Scilab program the implements the linear discrete heat exchanger simulation. The program below can be cut and pasted into Scilab’s editor.

// SIMPLE HEAT EXCHANGER CALCULATION FROM CONTROL GURU

// http://controlguru.com

// NEW SCILAB USERS SHOULD BUY

// Modeling and Simulation in Scilab/Scicos

// NOTE, EXPRESSIONS ENDING WITH SEMICOLONS WILL NOT PRINT WHEN EXECUTED

// EXPRESSIOBS WITHOUT SEMICOLONS WILL DISPLAY THE ASSIGNED VALUE(S)

Kp=-0.533; // PLANT GAIN, UNITLESS

Tp=1.3*60; // PLANT TIME CONSNTAT IN SECONDS

DT=round(0.8*60); // PLANT DEAD TIME IN SECONDS

T=1; // UPDATE TIME IN SECONDS

PVss=140-39*Kp; // STEADY STATE TEMPERATURE

// CALCULATE THE MODEL COEFFICIENTS

A=exp(-T/Tp); // DISCRETE TIME TRANSITION MATRIX WITH ONLY ONE ELEMENT

B=Kp*(1-A); // DISCRETE TIME INPUT COUPLING MATRIX

C=(1-A)*PVss; // OFFSET OR SYSTEM BIAS SO SYSTEM GOES TO PVss WITH NO CONTROL

// CALCULATE THE PI CONTROLLER GAINS USING THE IMC METHOD AS FOUND ON CONTROLGURU

Tc=max(Tp,8*DT); // DESIRED CLOSED LOOP TIME CONSTANT

Kc=Tp/(Kp*(Tc+DT)); // CONTROLLER PROPORTIONAL GAIN

Ti=Tp; // CONTROLLER INTEGRATOR TIME CONSTANT

for n=1:3600 // SIMULATE 60 MINUTES OF TIME. n = SECONDS

// INITIAL CONDITIONS, NOTE SCILAB ARRAYS START WITH 1

if n == 1 then

SP(1)=140; // INITIAL SET POINT

PV(1)=140; // PROCESS VARIABLE, TEMPERATURE = 140

CV(1)= 39; // CONTROL VARIABLE, OUTPUT % = 39%

CVi(1)=CV(1); // THE INTEGRATOR CONTRIBUTION TO THE OUTPUT. 39%

CVp(1)=0; // PROPORTIONAL CONTRIBBUTION TO THE OUTPUT

// ASSUME INITALY THERE IS NOT ERROR

end

// GENERATE THE SP AS A FUNCTION OF TIME

// IN THIS CASE THERE IS JUST A SIMPLE STEP CHANGE AS IN THE CONTROLGURU EXAMPLE

// ADD MORE SP CHANGES BY USING THE elsif KEYWORD

if n