fmincon question - does not converge to useful answer

Ok, so I am trying to write a program that can calculate the critical temperature and critical volume of a mixture of n components. This program is going to be part of a larger program that I am writing to do thermodynamic calculations.

This is how it works (or will work hopefully): The user inputs the required data on the components (critical temperature, critical pressure, mole fraction, and the component names). The program uses this data to minimize the Helmholtz free energy, thus calculating the critical point of the mixture. Basically, there are n+2 equations that equal zero to solve simultaneously (n number of PDE's for fugacity, 1 for the second PDE of the fugacity, and one constraining the perturbation vector to be always non-zero. I wrote a function that should do this, but I'm getting erroneous results. For the problem below, the critical temperature should be 462.64 K, the critical volume should be 228.83 cc/gmol, and the critical pressure should be 117 bar. The perturbation components should be 0.1137, 0.0485, and -0.9923, respectively. I'm having issues with getting fmincon converging to the right answer. If anyone can take a look at the code below and comment me, and my head that has been banging off the wall for two weeks would appreciate it.

The code that I wrote below only seems to vary the perturbation variables "del_n". It often returns the initial guesses for temperature and volume as the final answers. Is this the right way to approach this problem?

[Joshua]

function [Tc_mix,Pc_mix,Vc_mix,Zc_mix,exitflag,del_n] = EOS_crit(EOS,z,Pc,Tc,Zc,w,kij,method)
%This function will calculate the critical point parameters for a mixture
%of n components

%Origin paper can be found here:
%http://www.ijoticat.com/index.php/IJoT/article/viewFile/380/377

%Defined in another script, taken to be 83.14 cm^3.bar/mol.K
global R
%nc represents the number of components, this is defined as the length of
%one of the critical component row vectors (Tc, Pc, w); for example
%length(w)
global nc

%Feed composition normalization
%In the event of a typo, this ensures that the mol fractions sum to 1
z = z/sum(z(1:nc));

%Equation of state parameters based off of variable switch "EOS"
%This allows the user to select which EOS to use; this is defined in
%another script. Use 'PR' to evaluate
switch EOS
  case('vdW')
      sigma_1 = 0;
      sigma_2 = 0;
      Omega_a = 27/64;
      Omega_b = 1/8;
  case('RK')
      sigma_1 = 1;
      sigma_2 = 0;
      Omega_a = 0.42748;
      Omega_b = 0.08664;
  case('SRK')
      b0 = 0.480;
      b1 = 1.574;
      b2 = -0.176;
      sigma_1 = 1;
      sigma_2 = 0;
      Omega_a = 0.42748;
      Omega_b = 0.08664;
  case('PR')
      b0 = 0.37464;
      b1 = 1.54226;
      b2 = -0.26992;
      sigma_1 = 1 + sqrt(2);
      sigma_2 = 1 - sqrt(2);
      Omega_a = 0.457235;
      Omega_b = 0.077796;
end
%These are "scaling factors to scale the temp and volume to be closer in value to the mol fractions. Not currently used.
%VSF0 = Omega_b*R.*Tc./Pc;
%VSF = sum(z.*VSF0);
VSF = 1;
%TSF = 200;
TSF = 1;
Tc = Tc./TSF;
Pc = Pc./VSF;
b = Omega_b*R.*Tc./Pc;
b_mix  = sum(z.*b);

%Guesses for initial critical parameters
Tc_0 = sum(z.*Tc);
Vc_0 = (R*sum((z.*(Zc.*Tc)./Pc)));
del_n0 = z.^(2/3);

%del_n0 = [0.1137 0.0485 -0.9923];

%g=internal state vector of length nc+2
g0 = [Tc_0,Vc_0,del_n0];

%-------------------------------------------------------------------------%
%Switch fmincon calculation method based on user choice
switch method
  case(1)
      method = 'trust-region-reflective'; %Current default
  case(2)
      method = 'active-set';
  case(3)
      method = 'interior-point'; %(will become default in a future release)
  case(4)
      method = 'sqp';
end

%Set fmincon options
%options = optimset('Display','iter','MaxFunEvals',10000,'MaxIter',10000,'Algorithm',method,'TolX',1e-6,'TolCon',1e-6,'TolFun',1e-6);
options = optimset('FinDiffType','central','Display','iter-detailed','MaxFunEvals',10000,'MaxIter',10000,'Algorithm',method,'TolX',1e-6,'TolCon',1e-6,'TolFun',1e-6);
% lb = [];
% ub = [];
%Upper and lower constraints for iteration
%Tc_lb = lower limit for the mixture critical temperature
%Vc_lb = lower limit for the mixture critical volume
%Tc_ub = upper limit for the mixture critical temperature
%Vc_ub = upper limit for the mixture critical volume
%"del_n_lb/ub limit the vaules of del_n to [0,1] for component 1, and to
%[-1,1] for the other components
Tc_lb = min(Tc);
Vc_lb = (0.5*Vc_0);
Tc_ub = max(Tc);
Vc_ub = (2*Vc_0);
del_n_lb(1) = 0;
del_n_lb(2:nc) = -1;
del_n_ub(1) = 1;
del_n_ub(2:nc) = 1;
lb = [Tc_lb,Vc_lb,del_n_lb];
ub = [Tc_ub,Vc_ub,del_n_ub];

%Call fmincon to solve flash problem defined in "crit_param"
[g,exitflag] = fmincon('1',g0,[],[],[],[],lb,ub,@crit_param,options);

%Extract states for solution
Tc_mix = g(1)*TSF;
Vc_mix = g(2)*VSF;
del_n = g(3:end);   
Pc_mix = (R*Tc_mix/(Vc_mix - b_mix)) - (a_mix/((Vc_mix + sigma_2*b_mix)*(Vc_mix + sigma_1*b_mix)));
Zc_mix = (Pc_mix*Vc_mix)/(R*Tc_mix);

%-------------------------------------------------------------------------%
%Nested function to solve the energy calculation
function [c,ceq] = crit_param(g)

%Extract states:
%Critical temperature
Tc_mix = g(1);

%Critical volume
Vc_mix = g(2);

%Change in mol fraction
del_n(1:nc) = g(3:end);

%Physical constants
switch EOS
  case('vdW')
      a_Tr = ones(1,nc);
  case('RK')
      a_Tr = Tr.^(-1/2);
  otherwise
      m = b0 + b1.*w + b2.*w.^2;
      a_Tr = (1.0 + m.*(1.0 - sqrt(Tc_mix./Tc))).^2;
      for i=1:nc
          a(i) = (Omega_a*((R*Tc(i))^2)/Pc(i))*a_Tr(i);
      end
end

%Mixing rules
a_mix = 0; %b_mix = 0; 
aij = zeros(nc); a_a = 0;
for i=1:nc
  for j=1:nc
      aij(i,j) = sqrt(a(i)*a(j))*(1.0-kij(i,j));
      %"a_mix" mixture constant
      a_mix    = z(i)*z(j) * aij(i,j) + a_mix;
      a_a    = del_n(i)*del_n(j) * aij(i,j) + a_a;
  end
end
%-------------------------------------------------------------------------%

%Pre-allocate sums used below
alpha_sum = zeros(nc,nc);
part2_sum1 = zeros(nc);
part2_sum = zeros(1,nc);
alpha_k = zeros(1,nc);

%Constants required for the Helmholtz energy calculation
a_bar = a_a/a_mix;
K = Vc_mix/b_mix;
F1 = 1/(K-1);
F2 = (2/(sigma_1-sigma_2))*((sigma_1/(K+sigma_1))-(sigma_2/(K+sigma_2)));
F3 = (2/(sigma_1-sigma_2))*((sigma_1/(K+sigma_1))^2-(sigma_2/(K+sigma_2))^2);
F4 = (2/(sigma_1-sigma_2))*((sigma_1/(K+sigma_1))^3-(sigma_2/(K+sigma_2))^3);
F5 = (2/(sigma_1-sigma_2))*log(((K+sigma_1)/(K+sigma_2)));
F6 = F2-F5;
beta01 = b./b_mix;
beta01_bar = sum(del_n.*beta01);
N_bar = sum(del_n);

for i=1:nc
  for j=1:nc
      part2_sum1(i,j) = del_n(j)*aij(i,j);
      alpha_sum(i,j) = (z(i)*aij(j,i));
  end
  part2_sum(i) = sum(part2_sum1(i,:));
end
for i=1:nc
  alpha_k(i) = sum(alpha_sum(:,i))/a_mix;
end
alpha_bar = sum(del_n.*alpha_k);

%Nonlinear inequality constraints for fmincon
%c = Compute nonlinear inequalities at g - think about adding some!!!!

%f4 = Tc_mix-min(Tc);
%f4 = min(Pc)-(R*Tc_mix/(Vc_mix - b_mix)) - (a_mix/((Vc_mix + sigma_2*b_mix)*(Vc_mix + sigma_1*b_mix)));
c = []; 
%ceq = Compute nonlinear equalities at g
%Nonlinear equality constraints for fmincon (they equal zero)

%This represents the first order PDE of the component fugacity (i) with
%respect to component (j) at constant k,...nc * the pertubation changes
for l=i:nc
  for j=i:nc
      part2_sum1(i,j) = aij(i,j)*del_n(j);
  end
  dA(1,i) = ((R*Tc_mix)*(del_n(i)/z(i)+F1*(beta01(i)*N_bar+beta01_bar)+beta01(i)*F1^2*beta01_bar));...
            +((a_mix/b_mix)*((beta01(i)*beta01_bar*F3)-((F5/a_mix)*part2_sum(i))+(F6*(beta01(i)*beta01_bar-alpha_k(i)*beta01_bar-alpha_bar*beta01(i)))))
end

%This represents the second order PDE of the component fugacity (i) with
%respect to components (j) and (k) at constant m...nc
ddA = (R*Tc_mix)*(3*N_bar*(beta01_bar*F1)^2+3*(beta01_bar*F1)^3-sum((del_n.^3)./(z.^2)));...
    +(a_mix/b_mix)*(3*beta01_bar^2*(2*alpha_bar-beta01_bar)*(F2+F3+F6)-2*beta01_bar^3*F4-3*beta01_bar*a_bar*F6)

%Define nonlinerar equaliry constraints. All of the partial derivatives
%should equal zero at the critical point.
%f1 requires each if the individual derivatives to be equal to zero
%f2 requires that the first order PDE's sum to zero
%f3 requires that the second order PDE equal zero
%f4 requires that the del_n vector always be non-zero by normalizing it.
%f5 requires that the sum of the first order PDE's and the second order PDE
%be equal

f1(1:nc) = dA(1:nc)
f2 = sum(dA)
f3 = ddA
f4 = del_n*del_n'-1
f5 = sum(dA)-ddA
ceq = [f1';f2';f3';f4';f5']
end
end

%The critical data for the problem are:
%names = {'Methane','Nitrogen','n-Hexane'};
%Tc = [190.6 126.1 507.4];
%Pc = [46.04 33.94 30.12];
%w = [0.011 0.04 0.305];
%z = [0.46 0.07 0.47];
%mw = [16.043 28 86.177];
%Zc = [0.288 0.292 0.264];
%nc = length(w);
%kij(1,:) = [0 0.0278 0.0422];
%kij(2,:) = [0.0278 0 0.1496];
%kij(3,:) = [0.0422 0.1496 0];
%'method' = 3;

[Sean de Wolski]

What does the iter-detailed output look like?

[Joshua]

Optimization completed: The relative first-order optimality measure, 0.000000e+000,
is less than options.TolFun = 1.000000e-006, and the relative maximum constraint
violation, 7.976410e-007, is less than options.TolCon = 1.000000e-006.

Optimization Metric                                         Options
relative first-order optimality =  0.00e+000       TolFun =  1e-006 (selected)
relative max(constraint violation) =  7.98e-007    TolCon =  1e-006 (selected)

Mixture critical temperature, pressure calculated

exitflag = 1

First-order optimality measure was less than options. TolFun, and maximum constraint violation was less than options.TolCon

Results for mixture critical temperature and pressure:
The estimated critical pressure is 62.4587 Bar

The estimated critical temperature is 438.0623 K

The estimated critical volume is 376.2377 cc/mol

Change in molar fractions are:

Methane 0.95600632
Nitrogen -0.27234165
n-Hexane -0.10919044

As you can see, the method does converge given the constraints - it's just not converging to the answer I would expect.

[Joshua]

You know, I'm wondering if using fmincon is the way to go about solving this problem. I am looking for the values that result in the functions in the objective function equaling zero, not necessarily a minimum point. When the conditions are met, the objective function (the first and second partial derivatives of the fugacity calculation based on a cubic equation of state) should equal zero - within the tolerance level specified by the end user. It's possible that the minimum calculated by fmincon is a local minimum, and not a zero. Would it be better to use fsolve?

I have to admit, I'm not very good at using these functions. Most of the literature that covers this topic uses an iterative procedure based on a dampened Newton-Rhapson method to solve the system of PDE's. I figured why reinvent the wheel when Matlab has the tools already built in.

[Joshua]

I'm pretty sure that the initial values I am supplying are pretty good guesses. The value are based on molar fraction weighted averages - which are fairly good approximations for the hydrocarbon systems I am looking at (non-polar, no hydrogen bonding, ect). The fact that the function doesn't return the initial values when I use the actual answers as the initial guesses leads me to think that my methodology is wrong on this.

[Joshua]

 function [Tc_mix,Pc_mix,Vc_mix,Zc_mix,exitflag,del_n] = EOS_crit(input)
%This function will calculate the critical point parameters for a mixture
%of n components

%Origin paper can be found here:
%http://www.ijoticat.com/index.php/IJoT/article/viewFile/380/377

Tc = input(1,:);
Pc = input(2,:);
w = input(3,:);
z = input(4,:);
Zc = input(5,:);
kij = input(6:end,:);

global R
global nc
global EOS
global method

%nc represents the number of components, this is defined as the length of
%one of the critical component row vectors (Tc, Pc, w); for example
%length(w)

%Feed composition normalization
%In the event of a typo, this ensures that the mol fractions sum to 1
z = z/sum(z(1:nc));

%VSF0 = Omega_b*R.*Tc./Pc;
%VSF = sum(z.*VSF0);
VSF = 1;
%TSF = 200;
TSF = 1;
Tc = Tc./TSF;
Pc = Pc./VSF;

%Guesses for initial critical parameters
Tc_0 = sum(z.*Tc);
Vc_0 = (R*sum((z.*(Zc.*Tc)./Pc)));
del_n0 = z.^(2/3)/sum(z.^(2/3));

% Tc_0 = 462.64;
% Vc_0 = 228.83;
% del_n0 = [0.1137 0.0485 -0.9923];

%g=internal state vector of length nc+2
g0 = [Tc_0,Vc_0,del_n0]

%Switch fmincon calculation method based on user choice
switch method 
  case(1)
      method = 'trust-region-reflective'; %Current default
  case(2)
      method = 'active-set';
  case(3)
      method = 'interior-point'; %(will become default in a future release)
  case(4)
      method = 'sqp';
end

%-------------------------------------------------------------------------%

%Set fmincon options

%options = optimset('AlwaysHonorConstraints','bounds','FinDiffType','central','Display','iter-detailed','MaxFunEvals',10000,'MaxIter',10000,'Algorithm',method,'TolX',1e-12,'TolCon',1e-12,'TolFun',1e-12);
% lb = [];
% ub = [];
%Upper and lower constraints for iteration
%Tc_lb = lower limit for the mixture critical temperature
%Vc_lb = lower limit for the mixture critical volume
%Tc_ub = upper limit for the mixture critical temperature
%Vc_ub = upper limit for the mixture critical volume
%"del_n_lb/ub limit the vaules of del_n to [0,1] for component 1, and to
%[-1,1] for the other components
Tc_lb = min(Tc);
Vc_lb = (0.5*Vc_0);
Tc_ub = max(Tc);
Vc_ub = (1.2*Vc_0);
del_n_lb(1) = 0;
del_n_lb(2:nc) = -1;
del_n_ub(1) = 1;
del_n_ub(2:nc) = 1;
lb = [Tc_lb,Vc_lb,del_n_lb]
ub = [Tc_ub,Vc_ub,del_n_ub]
A = [];
b = [];
Aeq = [];
beq = [];

%Call fmincon to solve problem defined in "crit_func"
[g,exitflag] = fmincon(@crit_func,g0,A,b,Aeq,beq,lb,ub,@crit_param,options);

%Extract states for solution
Tc_mix = g(1)*TSF;
Vc_mix = g(2)*VSF;
del_n = g(3:end);   
Pc_mix = (R*Tc_mix/(Vc_mix - b_mix)) - (a_mix/((Vc_mix + sigma_2*b_mix)*(Vc_mix + sigma_1*b_mix)));
Zc_mix = (Pc_mix*Vc_mix)/(R*Tc_mix);
%-------------------------------------------------------------------------%

function [f,dA,ddA] = crit_func(g)
%This function will calculate the critical point parameters for a mixture
%of n components

Tc_mix = g(1);
Vc_mix = g(2);
del_n = g(3:end);

%Pre-allocate sums used below
alpha_sum = zeros(nc);
part2_sum1 = zeros(nc);
part2_sum = zeros(1,nc);
alpha_k = zeros(1,nc);
dA = zeros(1,nc);
aij = zeros(nc);

%Equation of state parameters based off of variable switch "EOS"
%This allows the user to select which EOS to use; this is defined in
%another script. Use 'PR' to evaluate
switch EOS
  case('vdW')
      sigma_1 = 0;
      sigma_2 = 0;
      Omega_a = 27/64;
      Omega_b = 1/8;
  case('RK')
      sigma_1 = 1;
      sigma_2 = 0;
      Omega_a = 0.42748;
      Omega_b = 0.08664;
  case('SRK')
      b0 = 0.480;
      b1 = 1.574;
      b2 = -0.176;
      sigma_1 = 1;
      sigma_2 = 0;
      Omega_a = 0.42748;
      Omega_b = 0.08664;
  case('PR')
      b0 = 0.37464;
      b1 = 1.54226;
      b2 = -0.26992;
      sigma_1 = 1 + sqrt(2);
      sigma_2 = 1 - sqrt(2);
      Omega_a = 0.457235;
      Omega_b = 0.077796;
end

%Physical constants
switch EOS
  case('vdW')
      a_Tr = ones(1,nc);
  case('RK')
      a_Tr = Tr.^(-1/2);
  otherwise
      m = b0 + b1.*w + b2.*w.^2;
      a_Tr = (1.0 + m.*(1.0 - sqrt(Tc_mix./Tc))).^2;
end

a = Omega_a.*a_Tr.*(((R*Tc).^2)./Pc);
b = Omega_b.*R.*Tc./Pc;

%Mixing rules
a_mix = 0; b_mix = 0; a_a = 0;
for i=1:nc
  for j=1:nc
      aij(i,j) = sqrt(a(i)*a(j))*(1.0-kij(i,j));
      a_mix    = z(i)*z(j) * aij(i,j) + a_mix;
      a_a    = del_n(i)*del_n(j) * aij(i,j) + a_a;
  end
  b_mix  = z(i)*b(i) + b_mix;
end

a_bar = a_a/a_mix;
K = Vc_mix/b_mix;
beta01 = b./b_mix;
beta01_bar = sum(del_n.*beta01);
N_bar = sum(del_n);

%Constants required for the Helmholtz energy calculation

F1 = 1/(K-1);
F2 = (2/(sigma_1-sigma_2))*((sigma_1/(K+sigma_1))-(sigma_2/(K+sigma_2)));
F3 = (2/(sigma_1-sigma_2))*((sigma_1/(K+sigma_1))^2-(sigma_2/(K+sigma_2))^2);
F4 = (2/(sigma_1-sigma_2))*((sigma_1/(K+sigma_1))^3-(sigma_2/(K+sigma_2))^3);
F5 = (2/(sigma_1-sigma_2))*log(((K+sigma_1)/(K+sigma_2)));
F6 = F2-F5;

for i=1:nc
  for j=1:nc
      part2_sum1(i,j) = del_n(j)*aij(i,j);
      alpha_sum(i,j) = (z(i)*aij(j,i));
  end
  part2_sum(i) = sum(part2_sum1(i,:));
end

for i=1:nc
     alpha_k(i) = sum(alpha_sum(:,i))/a_mix;
end

alpha_bar = sum(del_n.*alpha_k);

%This represents the first order PDE of the component fugacity (i) with
%respect to component (j) at constant k,...nc * the pertubation changes
for i=1:nc
  dA(i) = ((R*Tc_mix)*(del_n(i)/z(i)+F1*(beta01(i)*N_bar+beta01_bar)+beta01(i)*F1^2*beta01_bar))...
      +((a_mix/b_mix)*((beta01(i)*beta01_bar*F3)-((F5/a_mix)*part2_sum(i))+(F6*(beta01(i)*beta01_bar-alpha_k(i)*beta01_bar-alpha_bar*beta01(i)))));
end

%This represents the second order PDE of the component fugacity (i) with
%respect to components (j) and (k) at constant m...nc
ddA = (R*Tc_mix)*(3*N_bar*(beta01_bar*F1)^2+3*(beta01_bar*F1)^3-sum((del_n.^3)./(z.^2)))...
  +(a_mix/b_mix)*(3*beta01_bar^2*(2*alpha_bar-beta01_bar)*(F2+F3+F6)-2*beta01_bar^3*F4-3*beta01_bar*a_bar*F6);

dA_sum = sum(dA);

f = dA_sum - ddA;

end

function [c,ceq] = crit_param(g)
%This function will calculate the critical point parameters for a mixture
%of n components
% global nc
% global EOS
% Tc = g(3,:);
% Pc = g(4,:);
% w = g(5,:);
% z = g(6,:);
% kij = g(8:end,:);

[~,dA,ddA] = crit_func(g);

%Nonlinear inequality constraints for fmincon
%c = Compute nonlinear inequalities at g - think about adding some!!!!
c = []; 
%ceq = Compute nonlinear equalities at g  
%Define nonlinerar equaliry constraints. All of the partial derivatives
%should equal zero at the critical point.
%f1 requires each if the individual derivatives to be equal to zero
%f2 requires that the first order PDE's sum to zero
%f3 requires that the second order PDE equal zero
%f4 requires that the del_n vector always be non-zero by normalizing it.
%f5 requires that the sum of the first order PDE's and the second order PDE
%be equal

f1(1:nc) = dA(1:nc);
f2 = sum(dA);
f3 = ddA;
f4 = del_n*del_n'-1;
f5 = sum(dA)-ddA;
ceq = [f1';f2';f3';f4';f5'];

end
end

[Joshua]

For what it's worth, this is the code I've brewed up so far. I've nested the constraint function and the objective function in the main function.

[Changxu]

Hi, are you still there? I want to know the Tc, Pc and Vc calculated by PR EOS. Then know the deviation between critical properties obtained by PR EoS and the experimental data.

How can I achieve this? I don't know how to use Helmholtz energy calculationto realize it. And for now, I know how to calculate the saturated properties by using PR EoS.

Can you help me? Thanks.

Answer 1

[Matt J]

I don't think you've showed us your objective function. In your call to fmincon, you designate the objective function as '1'

fmincon('1',g0,[],[],[],[],lb,ub,@crit_param,options)

which should not work.

Aside from this, I'll point out that you have sqrt() operations all over the place. That will make the function non-differentiable at zero and undefined below 0, unless you're sure the quantites being sqrt-ed stay well above zero due to imposed bounds. Even with imposed bounds, you would probably need to make sure your bounds are always honored using the 'AlwaysHonorConstraints' option which is only available in certain methods (interior-point and sqp, I think).

[Joshua]

The objective function is nested in the main function. It's "@crit_param". The square roots are unfortunately unavoidable due to the nature of the equation of state used to calculate the pure species properties. They do stay above zero though. I'll give the 'AlwaysHonorConstraints' option a try and see if it makes any difference. From watching the code run, it looks like it does already do that - it just doesn't converge to the answer I think it should.

[Matt J]

No, that's not consistent with the code as you've shown it to us. crit_param is being used to specify constraints and is returning [c,ceq], not objective function values. You have a string '1' where the objective function should be.

Answer 2

[Matt J]

If you know what the solution should be, initialize FMINCON at that solution and see if FMINCON gives that solution back to you as its output. If it does, your original problem may have been a local minimum and you need to initialize better. If it does not, the solution you're expecting is just not a minimum of the function as you've provided it to fmincon. You may have an error translating your physics into code... a problem only you're equipped to address.

[Joshua]

If I plug the answers into the code as the starting points, the script still iterates on the variables "del_n", which are the perturbation values of each component. The final answer returns the same values as the initial guess, but returns different values for "del_n".

If the code was correct, there should be no iteration, since each function would be equal to zero. I've literally spent days looking at this code to see if there was an error, but I've yet to find one. I even modeled these equations in other software packages and compared the answers - they match up. I'm stumped.

[Matt J]

Forget about FMINCON then. Just concentrate on debugging your objective function. Plug in the unknown solution and run the objective function in the debugger. Step through line by line until you reach a line that doesn't give the expected result.