image thumbnail

Multivariate Data Analysis and Monitoring for the Process Industries

by

 

26 Feb 2008 (Updated )

Files from the webinar Multivariate Data Analysis and Monitoring for the Process Industries

ProcessMonitor
function ProcessMonitor

data = load('pcamodel');

numvars = size(data.varnames,2);

%% Create main figure
hfig = figure('units','normalized','position',[0.1,0.1,0.8,0.6],...
    'menubar','none','toolbar','none','nextplot','add',...
    'name','Process Monitor', 'NumberTitle', 'off',...
    'tag','MSPCmonitoringtool','DeleteFcn',@MSPCclosefcn);

    function MSPCclosefcn(hObject,eventdata) %#ok
        t = timerfind('Tag','readOPCtimer');
        stop(t);
        delete(t);
        da = getappdata(hObject, 'dahandle');
        disconnect(da);
        close(hObject);
    end

%% Store data in figure
setappdata(hfig, 'means', data.means);
setappdata(hfig, 'sds', data.sds);
setappdata(hfig, 'variances', data.variances);
setappdata(hfig, 'loadings', data.loadings);
setappdata(hfig, 'varnames', data.varnames);

%% Create 12 axes in the figure
for i = 1:12
    h(i) = axes('parent', hfig);
end

set(h,'units','normalized',...
    'XTick', [], 'YTick', [], 'Color', 'none',...
    'XColor', [0.5,0.5,0.5], 'YColor', [0.5,0.5,0.5],...
    'NextPlot', 'add', 'Box', 'on');

set(h,{'position'},...
    {[0.05,0.55,0.2,0.4];...
    [0.3,0.55,0.2,0.4];...
    [0.05,0.3,0.45,0.2];...
    [0.05,0.05,0.45,0.2];...
    [0.55,0.875,0.4,0.075];...
    [0.55,0.757,0.4,0.075];...
    [0.55,0.639,0.4,0.075];...
    [0.55,0.521,0.4,0.075];...
    [0.55,0.403,0.4,0.075];...
    [0.55,0.285,0.4,0.075];...
    [0.55,0.167,0.4,0.075];...
    [0.55,0.05,0.4,0.075]});

setappdata(hfig, 'axeshandles', h);

%% Add loadings to biplot
axes(h(1));
xlabel('Amount of Impurity','Color','k'); ylabel('Pattern of Impurity','Color','k'); xlim([-18,4]); ylim([-8,5]);
scatter(data.loadings(:,1)*data.scale(1),data.loadings(:,2)*data.scale(2),'ro','filled');
for i=1:numvars
    line([0,data.loadings(i,1)*data.scale(1)],...
        [0,data.loadings(i,2)*data.scale(2)],'Color','r');
    text(data.loadings(i,1)*data.scale(1)-2,data.loadings(i,2)*data.scale(2),data.varnames{i});
end

hdata(1) = scatter([],[],'b','filled');

%% Add control limits and labels to polar plot
axes(h(2))
xlim([-10,10]); ylim([-10,10]);
upperlimits = [4,3,-2,-3];
lowerlimits = [3,2,-3,-4];
colors = {'r','y','y','r'};
offset = 4;
for i = 1:4
    [X,Y]=pol2cart(...
        [0:1/99:1,1:-1/99:0]*2*pi,...
        [ones(1,100)* (upperlimits(i) + offset),...
        ones(1,100)*(lowerlimits(i) + offset)]);
    patch(X,Y,colors{i},'EdgeColor','none','FaceAlpha',0.4);
end
[X,Y]=pol2cart(...
    (0:1/numvars:1)*2*pi,...
    ones(1,numvars+1)*(4 + offset));
for i=1:numvars
    htext=text(X(i),Y(i),data.varnames{i});
    circlepos = (i-1)/numvars;
    if circlepos < 0.25
        set(htext,'VerticalAlignment','bottom','HorizontalAlignment','left');
    elseif circlepos == 0.25
        set(htext,'VerticalAlignment','bottom','HorizontalAlignment','center');
    elseif circlepos < 0.5
        set(htext,'VerticalAlignment','bottom','HorizontalAlignment','right');
    elseif circlepos == 0.5
        set(htext,'VerticalAlignment','middle','HorizontalAlignment','right');
    elseif circlepos < 0.75
        set(htext,'VerticalAlignment','top','HorizontalAlignment','right');
    elseif circlepos == 0.75
        set(htext,'VerticalAlignment','top','HorizontalAlignment','center');
    elseif circlepos < 1
        set(htext,'VerticalAlignment','top','HorizontalAlignment','left');
    end
end

hdata(2) = plot(0,0);

%% Add control limits to T2 and Residual plots

for i = [3,4]
    axes(h(i))
    xlim([0,730]);
    switch i
        case 3
            ylabel('T^{2}','Color','k');
            ylim([-5,100]);
            lims = data.tcrit;
        case 4
            xlabel('Observation Number','Color','k');
            ylabel('Residual','Color','k');
            ylim([0,6]);
            lims = data.distcrit;
    end
    upperlimits = [lims(1),lims(2)];
    lowerlimits = [lims(2),lims(3)];
    colors = {'r','y'};

    for j=1:2
        patch([0,0,730,730],...
            [upperlimits(j),lowerlimits(j),lowerlimits(j),upperlimits(j)],...
            colors{j},'EdgeColor','none','FaceAlpha',0.4);
    end
    
    hdata(i) = plot(0,0);
end


%% Add control lmits to variable plots
for i = 5:12
    axes(h(i));
    xlim([0,730]); ylim([-4,4]);
    if i==12
        xlabel('Observation Number','Color','k');
    end
    ylabel(data.varnames{i-4},'Color','k');
    upperlimits = [4,3,-2,-3];
    lowerlimits = [3,2,-3,-4];
    colors = {'r','y','y','r'};
    for j=1:4
        patch([0,0,730,730],...
            [upperlimits(j),lowerlimits(j),lowerlimits(j),upperlimits(j)],...
            colors{j},'EdgeColor','none','FaceAlpha',0.4);
    end
    
    hdata(i) = plot(0,0);
end

setappdata(hfig, 'plothandles', hdata);

%% Connect to OPC server
da = opcda('localhost', 'Iconics.SimulatorOPCDA.2');
connect(da);    
metalsgroup = addgroup(da);
additem(metalsgroup, 'Metals.Ag');
additem(metalsgroup, 'Metals.Ni');
additem(metalsgroup, 'Metals.Pb');
additem(metalsgroup, 'Metals.Bi');
additem(metalsgroup, 'Metals.Sb');
additem(metalsgroup, 'Metals.As');
additem(metalsgroup, 'Metals.Te');
additem(metalsgroup, 'Metals.Se');

%% Save handles to OPC connection and item group to figure
setappdata(hfig, 'dahandle', da);
setappdata(hfig, 'grouphandle',metalsgroup);

%% Create timer
t = timer('ExecutionMode','fixedRate','Period',0.5,'BusyMode','drop',...
        'Tag','readOPCtimer','TimerFcn', @readOPCtimerfcn);

%% Start timer
start(t);

%% Timer Callback
    function readOPCtimerfcn(hObject,eventdata) %#ok
        % Get handle to main figure
        hfig = findobj('Tag', 'MSPCmonitoringtool');
        % Get data and axes handles from figure
        h = getappdata(hfig, 'axeshandles');
        hdata = getappdata(hfig, 'plothandles');
        means = getappdata(hfig, 'means');
        sds = getappdata(hfig, 'sds');
        variances = getappdata(hfig, 'variances');
        loadings = getappdata(hfig, 'loadings');
        varnames = getappdata(hfig, 'varnames');
        metalsgroup = getappdata(hfig, 'grouphandle');
        numvars = size(varnames,2);
        % Read OPC data
        S = read(metalsgroup);
        metals = [S.Value];
        % Standardise data, and calculate PC scores, T2, and residuals
        metalstd = (metals-means)./sds;
        scores = metalstd*loadings(:,1:2);
        t2 = sum((scores(:,1:2).^2)./variances(1:2)',2);
        residuals = (metalstd - scores*loadings(:,1:2)');
        distance = sqrt(sum(residuals.^2,2)./(numvars-2));
        % Plot PC Scores on biplot
        xdata = get(hdata(1),'XData');
        ydata = get(hdata(1),'YData');
        colors = mat2cell([0,linspace(0.5,0.8,27);0,linspace(0.5,0.8,27);ones(1,28)]',ones(1,28),3);
        if size(xdata,2)<28
            xdata = [xdata, scores(:,1)];
            ydata = [ydata, scores(:,2)];
            colors = colors(1:size(xdata,2));
        else
            xdata = [xdata(2:end), scores(:,1)];
            ydata = [ydata(2:end), scores(:,2)];
        end
        set(hdata(1),'XData',xdata,'YData',ydata);
        patches = get(hdata(1),'Children');
        set(patches, {'MarkerFaceColor'},colors);
        % Plot variables on polar plot
        offset = 4;
        [X,Y]=pol2cart(...
            (0:1/numvars:1)*2*pi,...
            [metalstd(1,:),metalstd(1,1)] + offset);
        set(hdata(2),'XData',X,'YData',Y);
        % Plot T2
        xdata = get(hdata(3),'XData');
        ydata = get(hdata(3),'YData');
        xdata = [xdata,xdata(end)+1];
        ydata = [ydata,t2];
        set(hdata(3),'XData',xdata,'YData',ydata);
        % Plot Residual
        xdata = get(hdata(4),'XData');
        ydata = get(hdata(4),'YData');
        xdata = [xdata,xdata(end)+1];
        ydata = [ydata,distance];
        set(hdata(4),'XData',xdata,'YData',ydata);
        % Plot variables
        for k = 5:12
            xdata = get(hdata(k),'XData');
            ydata = get(hdata(k),'YData');
            xdata = [xdata,xdata(end)+1];
            ydata = [ydata,metalstd(k-4)];
            set(hdata(k),'XData',xdata,'YData',ydata);
            set(h(k),'XLim',[0 730],'YLim',[-4 4]);
        end
    end
        
end

Contact us