Support Articles
Test and Measurement
Automatic Step Testing with the OPC Toolbox
This demonstration shows an example of how to use the OPC Toolbox to automate step testing. A model of a methanol-water distillation column is run through a step test on the Reflux input to determine the open-loop system response of the top and bottom methanol concentrations.
The distillation column under test is used to separate a mix of methanol and water (the feed) into bottom products (mostly water) and a methanol-saturated distillate. The concentrations of methanol in the distillate and bottom products is denoted xD and xB respectively. Steam is used to vaporise the bottom liquids, and the vapour is mixed with the feed in the column. The vapour exiting the top of the column is condensed, and part of that condensate is fed back into the column as Reflux. By varying the steam flow rate and the reflux flow rate, concentrations xD and xB will vary.
In this example, we wish to determine a dynamic model of the response of the column to variations in the Reflux flow rate, by fitting a first order model to the responses of the outputs to changes in the Reflux input.
We will perform this task completely automatically, using the OPC Toolbox.
The "actual" plant used in this example is a Simulink model of the distillation column. The plant is modeled as a set of first order responses with time delays. Rate limits are induced on the inputs to the column.
The model of the column is adapted from Wood, R.K. and M.W. Berry, "Terminal Composition Control of a Binary Distillation Column," Chemical Engineering Science, 28 (1973), pp. 1707-1717, and is the same model used in the Time Delays example in the Control System Toolbox.
This example utilises a specially configured Simulation Server from ICONICS, to provide the tags for the distillation column. For step testing, we define some tags to log outputs, and to write input step signals. The tags are defined in variables as we use them frequently in the step test.
da = opcda('localhost', 'ICONICS.Simulator');
connect(da);
tagReflux = 'Distillation.Reflux';
tagSteam = 'Distillation.Steam';
tagXD = 'Distillation.xD';
tagXB = 'Distillation.xB';
Since we will be stepping the reflux value in this test, we create a group containing that MV tag.
grpWrite = addgroup(da, 'WriteGroup');
itmWrite = additem(grpWrite, tagReflux, 'double');
The PVs are the top and bottom distillates. These are represented on our OPC server by the tags xD and xB respectively. In order to log data from these tags, we create a logging group. We must include the MV in that logging group.
grpLogging = addgroup(da, 'LoggingGroup');
itmLogging = additem(grpLogging, {tagXD, tagXB, tagReflux}, repmat({'double'},1,3));
This logging task must log data for as long as the step test needs to be performed. A Records Acquired callback function is used to drive the MV with the desired signal, at the required time.
In this particular case, we will use a zero-mean double-pulse signal, with the second pulse three times longer than the first. This excites sufficient modes in the system without causing too many process changes. To make up that signal, we use the
sampleTime = 0.5; % Half-second sampling
pulseWidth = 20; % A 20-second pulse width
pulseDelta = 0.1; % A 10% change in the reflux level
pulseRatio = 3; % Duration of long pulse relative to short one.
We can now derive the pulse train, defined using the times in the previous section.
pulseTime = cumsum([0, 1, 1, 1, pulseRatio, pulseRatio, pulseRatio]*pulseWidth);
pulseData = [0, 1, -1, 1, -1, 0, 0]*pulseDelta;
stairs(pulseTime, pulseData);ylim([-1, 1]*pulseDelta*1.1);
title('Test signal for exciting Reflux input.');

Given the specification of the test pulse, we can determine the logging duration, and configure the function that will be used to drive the inputs to the process with that pulse.
numSamples = ceil(pulseTime(end)./sampleTime);
set(grpLogging, 'UpdateRate', sampleTime, 'RecordsToAcquire', numSamples);
set(grpLogging, ...
'RecordsAcquiredFcn', {@autostepMVDriver, grpWrite, pulseTime./sampleTime, pulseData}, ...
'RecordsAcquiredFcnCount', floor(pulseWidth/sampleTime));
The RecordsAcquiredFcn properties defines a callback (function) to execute after every RecordsAcquiredFcnCount records have been logged. We use the function to write the step test pulse contained in pulseData to the distillation column's Reflux input. The function autostepMVDriver code is shown below.
type autostepMVDriver
function autostepMVDriver(grp, event, grpWrite, pulseSample, pulseData)
%AUTOSTEPMVDRIVER Writes pulse train for AutoStep test
ind = (pulseSample==event.Data.RecordsAcquired);
if any(ind)
writeasync(grpWrite, pulseData(ind));
end
disp(sprintf('Logging Data: %d/%d records.', event.Data.RecordsAcquired, grp.RecordsToAcquire));
This step ensures that the process is at a reasonable steady-state prior to commencing with the step test.
write(grpWrite, pulseData(1));
pause(15);
This next few lines of code start the logging process, then wait for that process to complete.
disp('Logging starts...');
start(grpLogging)
wait(grpLogging)
Logging starts...
Logging Data: 40/480 records.
OPC WriteAsync event occurred at local time 21:42:24
Transaction ID: 7
Group Name: WriteGroup
1 items written.
Logging Data: 80/480 records.
OPC WriteAsync event occurred at local time 21:42:44
Transaction ID: 8
Group Name: WriteGroup
1 items written.
Logging Data: 120/480 records.
OPC WriteAsync event occurred at local time 21:43:04
Transaction ID: 9
Group Name: WriteGroup
1 items written.
Logging Data: 160/480 records.
Logging Data: 200/480 records.
Logging Data: 240/480 records.
OPC WriteAsync event occurred at local time 21:44:04
Transaction ID: 10
Group Name: WriteGroup
1 items written.
Logging Data: 280/480 records.
Logging Data: 320/480 records.
Logging Data: 360/480 records.
OPC WriteAsync event occurred at local time 21:45:04
Transaction ID: 11
Group Name: WriteGroup
1 items written.
Logging Data: 400/480 records.
Logging Data: 440/480 records.
Logging Data: 480/480 records.
OPC WriteAsync event occurred at local time 21:46:04
Transaction ID: 12
Group Name: WriteGroup
1 items written.
The OPC engine logged the distillation column input (Reflux) and outputs (top and bottom concentrations) during the step test. To retrieve the data, we specify that the values must be returned as a double matrix. With the values, we get the item IDs, qualities and time stamps for each value.
[itmID, val, qual, tStamp] = getdata(grpLogging, 'double');
The OPC Specification does not guarantee that items returned from the server are in the same order as the client requested them. Therefore, we have to find the indices of the items in the returned data.
indMV = strmatch(tagReflux, itmID);
indXD = strmatch(tagXD, itmID);
indXB = strmatch(tagXB, itmID);
We will plot all the data in separate sub-plots.
timeVec=tStamp(:,indXD);
subplot(3,1,3)
stairs(timeVec, val(:,indMV));
ylabel('MV (Reflux)'); xlabel('Time (Absolute)');
datetick x keepticks
xL = xlim;
subplot(3,1,1)
stairs(timeVec, val(:,indXD));
ylabel('xD');
set(gca,'XTick',[]);
xlim(xL);
title(sprintf('Automated Step Test started %s', datestr(tStamp(1,indMV))));
subplot(3,1,2)
stairs(timeVec, val(:,indXB));
ylabel('xB');
set(gca,'XTick',[]);
xlim(xL);

Since we are finished with the step testing task, we can disconnect the client from the OPC Server and delete the OPC Toolbox objects. Deleting the client deletes all of the client's children also.
disconnect(da);
delete(da);
To perform system identification, we must have the System Identification Toolbox installed.
if ~license('test','Identification_Toolbox')
disp('The rest of this demo requires the System Identification Toolbox.');
return;
end
Begin by creating data sets. The SysID Toolbox cannot model multi-output process models, so we create two data sets.
idXD = iddata(val(:,indXD)-val(1,indXD), val(:,indMV), sampleTime, ...
'OutputName', 'xD', 'InputName', 'Reflux');
idXB = iddata(val(:,indXB)-val(1,indXB), val(:,indMV), sampleTime, ...
'OutputName', 'xB', 'InputName', 'Reflux');
Now fit a delayed first order and a delayed second order underdamped model to the xD (top) output and compare the results.
prXD1 = pem(idXD, 'P1D', 'InitialState', 'zero', 'Name','prXD1');
prXD2 = pem(idXD, 'P2D', 'InitialState', 'zero', 'Name','prXD2');
figure;
compare(idXD, prXD1, prXD2, 'InitialState', 'zero');
Warning: Td.max has been set to 15 sec (30 samples).
Warning: Td.max has been set to 15 sec (30 samples).

From the results, the first order model has a reasonable fit. However, the second-order model seems to capture the faster dynamics better.
present(prXD2);
Process model with transfer function
K
G(s) = ------------------ * exp(-Td*s)
(1+Tp1*s)(1+Tp2*s)
with K = 12.199+-0.48149
Tp1 = 13.692+-2.3778
Tp2 = 5.9076+-1.8069
Td = 3.3895+-0.57923
Estimated using PEM from data set idXD
Loss function 0.00353146 and FPE 0.00359081
Created: 21-Sep-2004 21:46:20
Last modified: 21-Sep-2004 21:46:23
For comparison, the distillation model used in this example is a linear approximation, with the reflux to xD path modelled as a first-order system with a time constant of 16.7 seconds, a gain of 12.78, and a time delay of 1 second. The model mismatch is due to the addition of a rate limiter in the reflux input stream, which introduces additional dynamics.