Attempting to make a plot with a cut off.

I am attempting to make a graph looking like the attached screen shot using the data I have collected. The data is stored in Costs in the code below I juat don't know how I can create this graph:

function CostOptimization()

   %Define constant coefficients

   C1 = 0.04;

   C2 = 0.11;

   C3 = 0.17;

   psat1 = 0.193;

   psat2 = 0.085;

   psat3 = 0.061;

   % Constraint barring for each cost and component (Constraints from
   % pressure, and minimum composition)

   Pboundx2 = [];

   PboundC = [];

   xboundx2 = [];

   xboundC = [];

   Pxboundx2 = [];

   PxboundP = [];

   % Calc all possible costs from all combinations of x_1, x_2, and x_3

   % At a resolution of 0.001 steps in mole fractions

   Costs = zeros(701,701,5);

   for i = 1:1:701

   for j = 1:1:701

   % Calculate mole fractions

   x_1 = 0.1 + (i-1)*0.001;

   x_2 = 0.1 + (j-1)*0.001;

   x_3 = 1 - x_1 - x_2;

   % Do not calc values for when the fractions do not sum to 1

   if (x_3 < 0.0999)

   Costs(i,j,:) = NaN;

   continue

   end

   % Calculate values of activity coeffiecents at current point

   [gamma_1, gamma_2, gamma_3] = CalcGammas(x_1, x_2, x_3);

   % Calculate vapor pressure by combining vapor pressures activity
   % coefficents and compositions of each component

   P = gamma_1 * x_1 * psat1 + gamma_2 * x_2 * psat2 + gamma_3 * x_3 * psat3;

   % Calculate cost based on given equation

   C = C1 * x_1 + C2 * x_2 + C3 * x_3 + (2.5 * P.^2);

   % Store the cost and vapor pressure at the current mole fracs

   Costs(i,j,1) = x_1;

   Costs(i,j,2) = x_2;

   Costs(i,j,3) = x_3;

   Costs(i,j,4) = P;

   Costs(i,j,5) = C;

   % Store points that are along the x_3 = 1 boundary

   if (x_1+x_2==.9) && P<=.15

   xboundx2 = [xboundx2 x_2];

   xboundC = [xboundC C];

   end

   if (x_1+x_2==.9)

   Pxboundx2 = [Pxboundx2 x_2];

   PxboundP = [PxboundP P];

   end

   % Store points along the boundary where pressure exceeds 0.15

   if abs(P-.15) < 0.00001

   Pboundx2 = [Pboundx2 x_2];

   PboundC = [PboundC C];

   end

   % Do not store cost if pressure is too high

   if P > 0.15

   Costs(i,j,5) = NaN;

   end

   end

   end

   % Find minimum cost in cost matrix

   min_Cost = min(min(Costs(:,:,5)))

   for i = 1:1:701

   for j = 1:1:701

   if Costs(i,j,5) == min_Cost

   min_x_1 = Costs(i,j,1)

   min_x_2 = Costs(i,j,2)

   min_x_3 = Costs(i,j,3)

   break

   end

   end

   end

end

%

Answer 1

[Walter Roberson]

Never compare for floating point equality between two values.

Exception: if it is known that one of the values is an exact copy of what the target value is. For example it is okay to use

x = rand(1,50);
mx = min(x);   %extract something from the array
same_as_min = (x == mx);

because in this csse the value to be compared is promised to be bit-for-bit identical.

In your code you have tests of the form

x1+x2==0.9

Immediately you can see that the calculation (x1+x2) is not extracting a bit-for-bit copy of a value: it is calculating a new value. And then you are trying to compare the new value to whatever MATLAB happens to decide to represent 0.9 as this time around.

If you think of every floating point constant such as 0.9 that you wrote in code as producing a random value close to what you want, and think of it producing a different random value every time you write the constant then you will get closer to what happens inside the programming language. This applies for every floating point constant that does not have integral value, except that the floating point constants that are an integral value plus 1/2, 1/4, 3/4, 1/8, 3/8, 5/8, 7/8 are safe. (Some others too.) But 0.1, 0.2, 0.3, 0.4, 0.6, 0.7, 0.8 and 0.9 are not safe.

Each new calculation of x1+x2 should be treated as if it produces a random value "close to" the "real" value.

So what should you do? You should be calculating whether the sum is "close enough" to your target value. For example

if abs(x1 + x2 - 0.9) < 10^(-6)

This is exactly the reason that your code already has

if abs(P-.15) < 0.00001

rather than the (not good) "if P==0.15" -- the code is recognizing that P only needs to be "close enough" to being 0.15. If it were to try to test for equality, it would probably not match exactly due to round-off problems.