%
% Function to calculate operating costs about a
% constraint assuming Normal distributions
%
%     [ h, cost, sf ] = costcon( oldmean, oldstd, 
%               redn, coef, constr [, 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
%         constr  = process variable constraint value
%         optim   = optimal value if not oldmean
%
% Bob Newell, May 1996
%
%----------------------------------------------------
%
function [ h, cost, sf ] = costcon( oldmean, ...
                 oldstd, redn, coef, constr, optim )
%
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;
m = constr - redn * ( constr - m );
v = redn * oldstd * redn * oldstd;
newd = exp( - (x-m) .* (x-m) / 2 / v ) ...
                          / sqrt( 2 * pi * v );
plot( x, oldd, '--', x, newd, '--' );
ax = axis;
line( [ constr constr ], [ 0 ax(4)/2 ] );
xlabel( 'Process Variable' );
ylabel( 'Scaled Values' );
hold on
% generate cost curves
c = polyval( coef, abs(x-constr) );
sf(1) = ax(4) / max(c);
plot( x, c * sf(1), ':' );
% 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), '-' ); 
% calculate cost curve areas
if constr > oldmean,
  pos = oldmean - 3 * oldstd;
else
  pos = oldmean + oldstd;
end;
cost(1) = sum( oldoc ) * ( b - a ) / 100;
text( pos, 0.9 * ax(4), ...
               [ 'Old Cost = ' num2str( cost(1) ) ] );
cost(2) = sum( newoc ) * ( b - a ) / 100;
text( pos, 0.8 * ax(4), ...
               [ 'New Cost = ' num2str( cost(2) ) ] );
%
hold off
return;
%------------------------------------------------------
% the end
