Using fmincon to fit experimental data to simulated data

Previously I have been using fminsearch to fit my simulated data to experimental data for my PhD project. The results have been good up until now.

"The code is at the bottom of this post"

I have found a new dataset which requires a constraint to be added to my previous equation. In other words, to ensure that my model behaves as necessary, the parameters must not fall below some threshold value.

I have N = Nmax (1 - ((E1 + C)/(C + E2))).

Due to an initial condition, the above equation must not drop below the aforementioned threshold value. Therefore, the second term (fraction one) in the equation above must not fall below 0.938 (random value).

Thus, ((E1 + C)/(C + E2)) <= 0.938

Rearranging this, I can derive an expression for E1 in terms of E2.

Due to the constraints on the parameters here, fminsearch surely would not be applicable? I therefore tried to use fmincon, and meeting the syntax requirements, I had A = [C -0.938] and b = 0.938*C (obtained from rearranging the above inequality - E1*C - 0.938E2 <= 0.938*C).

In a previous post the fit was rather weird, but I overcome that issue. Now, the fit begins and manages to go through several iterations acheiving fairly good fits throughout. However, after sometime, the solver seems to fail (I assume this is because of the parameters selected by fmincon) and the code fails before its complete.

CODE:

%Model:
%dc/dt = gc(1-c/K); where K = Kmax*(1-(Kmax*cd)/(Kcd + cd))

function FittingDrug_Wis_Model12v2_060220

clc
clear all
clear workspace

%Initial guesses
Kmax = 0.964; 
Kcd = 1.12e-8; 

p0 = [Kmax Kcd];
%Constraints for fmincon
A = [4e-8 -0.938];
b = 0.938*4e-8;


%LSQCURVE
 %Initial guess vector
% options = optimoptions('fmincon','Display','iter','Algorithm','sqp');
% x = fmincon(@(par)ExpData(par),p0,A,b,options)
% x = fmincon(@(par)ExpData(par),p0,Aeq,beq,options)
x = fmincon(@(par)ExpData(abs(par)),p0,A, b)

function err=ExpData(param)
%% Data from paper and other sources
format long %Outputted values have many dp's
Kmax = param(1); %Extract V from params vector
Kcd = param(2); %Extract Pibar from params vector

%Known values:
g = 0.55072;
K = 1.162e6;
Cd = 4e-8;

P = [g K Cd];

%Dimensional time in days
texpw = [0 5 8 10 13 15 20];

%Med drug data
cm = [7.17e4 2.29e5 2.65e5 2.68e5 2.97e5 2.95e5 2.37e5];
errm = [NaN NaN NaN NaN NaN NaN NaN];

%% Solver set-up and data extraction
%Initially, we will only be fitting for a single (high) drug
%initial values for c and pi
ch0 = cm(1); 

%Solver for fitting (high drug)
[t,xh] = ode15s(@(t,x)odeDrug(t,x,P),texpw,ch0);

%Data from fitting
cellND = xh(:,1); %Extracting the cell data from the solver * by K
cellND = transpose(cellND);

%% Error calculation
%calculate difference between simulated and experimental value
err=0; %Counter set-up
for j=1:length(texpw)
    err = err + ((cellND(j) - cm(j))/cm(j))^2;
end

%% Plotting

%'Live' plot of fitting process
figure(2)
plot(t,cellND, 'k', 'LineWidth',2.5)
hold on
scatter(texpw, cm, 'filled', 'g', 'Linewidth',2.5)
errorbar(texpw, cm, errm, 'g', 'LineStyle','none','Linewidth',1.5);
xlabel('Time (Days)'), ylabel('Cell Number'), title('Fitting for the high drug')
set(gca,'fontsize',32)
drawnow

%% Solvers

function dxdt = odeDrug(t,x,P)
%The coupled non dimensional equations will do in here.

g = P(1);
K = P(2);
Cd = P(3);

dxdt = x(1)*g*(1 - (x(1)/(K*(1 - ((Kmax*Cd)/(Cd+Kcd))))));
   
end
end
end