%
% Function to calculate operating costs about an
% optimum assuming Normal distributions
%
%     [ h, cost, sf ] = costopt( oldmean, oldstd, 
%               redn, coef [, optim ] )
%
% where   h       = plot handle
%         cost    = vector of old and new costs
%         sf      = vector of plot scaling factors
%                   for cost and operating cost
%         oldmean = mean of old performance
%         oldstd  = std of old performance
%         redn    = std reduction factor
%         coef    = vector polynomial of cost of
%                   deviating from the optimum
%                   cost = coef(n) + coef(n-1)*devn +
%                           coef(n-2)*devn^2 + ...
%                   where devn is deviation of process
%                   variable from the optimum
%         optim   = optimal value if not oldmean
%
% Bob Newell, May 1996
% version for control book, Feb 1997
%
%----------------------------------------------------
%
function [ h, cost, sf ] = costopt( oldmean, ...
                 oldstd, redn, coef, optim )
%
fntsz = 24;
h = figure;
% determine the x axis range
b = oldmean + 4 * oldstd;
a = oldmean - 4 * oldstd;
x = a:(b-a)/100:b;
% generate Normal distributions
m = oldmean;
v = oldstd * oldstd;
oldd = exp( - (x-m) .* (x-m) / 2 / v ) ...
                          / sqrt( 2 * pi * v );
if nargin > 4,
  m = optim;
end;
v = redn * oldstd * redn * oldstd;
newd = exp( - (x-m) .* (x-m) / 2 / v ) ...
                          / sqrt( 2 * pi * v );
plot( x, oldd, '-', x, newd, '-' );
set( gca, 'fontsize', fntsz );
xlabel( 'Process Variable' );
ylabel( 'Distributions and Cost Penalty' );
hold on
% generate cost curves
c = polyval( coef, abs(x-m) );
ax = axis;
sf(1) = ax(4) / max(c);
plot( x, c * sf(1), '--' );
hold off
%
%
h = [ h figure ];
% calculate opcost curves
oldoc = oldd .* c;
newoc = newd .* c;
sf(2) = ax(4) / 4 / max( [ max(oldoc) max(newoc) ] );
plot( x, oldoc * sf(2), '-', x, newoc * sf(2), '-' ); 
set( gca, 'fontsize', fntsz );
xlabel( 'Process Variable' );
ylabel( 'Operating Costs' );
% calculate cost curve areas
cost(1) = sum( oldoc ) * ( b - a ) / 100;
% text( m + oldstd, 0.9 * ax(4), ...
%               [ 'Old Cost = ' num2str( cost(1) ) ] );
cost(2) = sum( newoc ) * ( b - a ) / 100;
% text( m + oldstd, 0.8 * ax(4), ...
%               [ 'New Cost = ' num2str( cost(2) ) ] );
%
hold off
return;
%------------------------------------------------------
% the end
