fmincon stopping/taking too long after initial point

Hello, I'm trying to find optimal parameters for a cancer cell division model that involves ode's and some constraints ('feedbacks').

I'm trying to fit a set of data (cell growth over 20 days, measured each day although some points are missing) to the curve. I have to use fmincon, with the objective function being the sum of squared error, and not lsqcurve because of the nonlinear constraints included.

Here's my code, excluding the bounds and linear constraints since those are just numbers:

function pfit
%%x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)
% k: p1 / p2 / q1 / q2 / b1 / b2 / g1 / g2 / g3 / g4 / v1 / v2 / d1 / d2 / d3
%%setup
% data
td = 0:1:20; %time
td = td(:);
xd = [0.15, 0.1, 0.11, 0.16, 0.18, 0.16, 0.4, 0.5, 1, 1.5, NaN, NaN, 2.9, 3.1, 3, NaN, NaN, NaN, NaN, NaN, NaN]* 10^6; % data over 20 days
k0 = [0.5, 0.5, 0.2, 0.1, 8e-11,4e-12, 10e-12, 10e-12, 10e-13, 10e-12, 0.5, 1, 5e-4, 0.0025, 0.005];
options = optimoptions('fmincon','Algorithm','sqp','Display', 'iter','FiniteDifferenceStepSize',1e-3,'StepTolerance',1e-5);

%%nonlinear inequality(c) and equality (ceq) constraint functions
nonlcon = @nlcon;
    function [c,ceq] = nlcon(k)
        ceq = (k(1)-k(3))*k(11)-k(13);
        c = -k(11)*(k(1)-k(3))+k(14);
    end
%%ode setup
    function dx = tumor (t,x,k,tau)
        k(11)= k(11)/(1+k(5)*(x(3)*(t-tau))^2); %type 1 feedback
        k(12) = k(12)/(1+k(6)*(x(3)*(t-tau))^2);

          k(1) = k(1)/(1+k(7)*((x(3)*(t-tau))^2)); %type 2 feedback
          k(3) = k(3)/(1+k(8)*((x(3)*(t-tau))^2)); %type 2 feedback
          k(2) = k(2)/(1+k(9)*((x(3)*(t-tau))^2)); %type 2 feedback
          k(4) = k(4)/(1+k(10)*((x(3)*(t-tau))^2)); %type 2 feedback

          dx = zeros (3,1); %x(1)= CSC, x(2) = PC, x(3) = TDC
          dx(1) = (k(1)-k(3))*k(11)*x(1)-k(13)*x(1);
          dx(2) = (1-k(1)+k(3))*k(11)*x(1)+(k(2)-k(4))*k(12)*x(2)-k(14)*x(2);
          dx(3) = (1-k(2)+k(4))*k(12)*x(2)-k(15)*x(3);
      end
  %%solve ode
      function err = tumorode(k) %least square error
          xinit = [10^5,10^5,10^5];
          tau = 2;
          tspan = [0:0.001:20];

          [t,x]= ode45(@(t,x) tumor(t,x,k,tau),tspan,xinit);
          x = x(:,1) + x(:,2) + x(:,3);
          x_fit = [];

          for i = 1:21
              x_fit = [x_fit, x(t == td(i))]; %match ode45 data with actual data for each day 
          end
          x_fit = x_fit(~isnan(x_fit));
          xd = xd(~isnan(xd));
          err = sum(sum((x_fit(:) - xd).^2));
      end

%%solve
[k,fval] = fmincon(@tumorode, k0, A, b, [], [], lb, ub, nonlcon, options);
end

When I run fmincon, I just get

   Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                         step    optimality
      0          16    3.798016e+14     1.495e-01     1.000e+00     0.000e+00     7.661e+16

and from there, matlab is 'busy' for an indefinite amount of time without ever moving on. I have loosened the step tolerance/size but it doesn't seem to be working (or perhaps I have done it incorrectly). Could you tell me what I can do (rearrange the code, fix options, etc) to avoid this lag? I have a feeling that the ode45 setup / loops are slowing things down, but I am a newcomer to this program/toolbox and don't know how to debug this correctly.

Any help or advice on this specific problem or optimization w/ fmincon in general would be appreciated. (+ I only have access to this toolbox)

-SY

Answer 1

[Matt J]

We can't run your code, but I suspect it's because your ODE has a singularity at t=tau and your code tries to integrate right through it, as far as I can tell. This will generate lots of NaNs and NaNs make things really slow.

[SY]

Hi Matt, thank you for the answer. You are right that I'm running into a singularity at t = tau. I tried using ode15s instead of ode45 to jump over this point, and the code now does execute more than one iteration. However, fmincon is halting prematurely:

 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0          16    3.794696e+14     1.495e-01     1.000e+00     0.000e+00     4.862e+21  
Warning: Failure at t=2.362015e-11.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (8.391562e-26) at time t. 
> In ode23s (line 401)
  In pfit/tumorode (line 75)
  In sqpInterface
  In fmincon (line 810)
  In pfit (line 88) 
    1         125    3.792989e+14     1.495e-01     3.928e-15     2.518e-01     5.550e+21  
    2         141    5.077972e+14     2.289e-03     1.000e+00     1.065e-01     1.451e+22  
    3         157    3.910419e+14     2.238e-05     1.000e+00     9.866e-03     5.277e+21  
    4         173    3.910416e+14     3.681e-11     1.000e+00     1.649e-05     1.116e+17  
    5         189    3.910416e+14     1.388e-17     1.000e+00     2.713e-11     9.101e+14  
Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

Clearly fmincon is stopping at a point where fval is not minimized. How can I troubleshoot this?

Thanks,

SY

[Matt J]

Clearly fmincon is stopping at a point where fval is not minimized.

How is that clear?

[SY]

From what I understand, the third column of the iterative display indicates the value of the objective function --- 4e14 doesn't seem to be the minimum value possible. Is that right?

[Matt J]

Why not? The first-order optimality dropped 7 orders of magnitude from the initial point to the final iteration. What does the exitflag say?

[x,fval,exitflag,output]=fmincon(...)

[SY]

The exitflag = 1, but when I plug the parameters that the code gives me at the final iteration into my equation, the results don't match my data points very well. I'm suspecting that there is a problem with a part of my code other than fmincon (which I now understand better). Thank you for your help on my question, I really appreciate it.